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CHAPTER  1:  INTRODUCTION 


BACKGROUND 

Airfield  pavement  design  is  a  complex  blend  of  relatively  simple  linear  elastic  theory, 
fatigue  concepts,  correlations  with  small  and  full-scale  tests,  and  pragmatic  adjustments  to 
reflect  observations  of  in-service  pavements.  This  philosophy  served  the  design  community 
well  for  many  years  as  it  allowed  total  thickness,  asphalt  concrete  pavement  thickness,  and 
material  requirements  for  constituent  layers  in  the  pavement  to  be  determined  to  avoid  a  pre¬ 
selected  level  of  distress  in  the  pavement.  For  airfields,  this  level  of  distress  at  “design” 
failure  was  selected  to  be  one  inch  of  shear  rutting  in  the  subgrade  or  fatigue  cracking  of  the 
asphalt  concrete. 

However,  today’s  designers  are  being  asked  to  predict  pavement  performance.  This 
is  a  far  more  complex  task  than  simply  providing  safe  thickness  and  specifications  for  the 
material.  To  deal  with  this  new  challenge,  the  design  community  must  have  material  models 
that  predict  cumulative  deformations  under  repetitive  aircraft  loads.  With  heavy  loading, 
such  as  may  be  encountered  with  many  airfields,  the  nonlinear  response  of  base  course 
materials  must  be  considered  when  predicting  pavement  performance.  The  advances  made  in 
computational  mechanics  have  created  new  tools  of  application  for  this  type  of  problem. 
Theoretically  rigorous  material  models  may  be  implemented  within  many  of  the  general- 
purpose  finite  element  computer  programs  available  today.  In  order  to  apply  these  material 
models,  mechanical  response  data  is  required  to  calibrate  the  necessary  model  parameters 
(Barker  and  Gonzalez,  1991). 

OBJECTIVE 

The  objective  of  this  research  was  to  provide  an  analytical  method  for  modeling  the 
response  of  unbound  granular  layers  in  flexible  pavements  subjected  to  aircraft  loads.  The 
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essential  features  of  pavement  response  that  are  required  from  a  constitutive  model  include 
non-linear  elastic  response,  permanent  or  plastic  deformation  after  yield,  cyclic  loading,  strain 
softening/hardening,  and  shear  dilatancy.  A  pavement  model  should  be  simple  in  operation, 
implementation  and  calibration.  The  model  must  be  executable  within  a  proven  general 
purpose  finite  element  code  like  ABAQUS  from  HKS,  Inc.  The  model  must  also  provide 
pavement  analysts  with  the  capability  of  predicting  the  performance  of  unbound  materials 
under  traffic  loadings. 

ORIGINALITY 

The  contribution  or  originality  of  the  research  is  in  the  following  area:  The 
identification,  implementation,  and  evaluation  of  a  new  constitutive  model  that  can  provide  for 
response  predictions  of  stresses  in  granular  pavement  layers  for  current  and  future  aircraft. 

SCOPE  OF  WORK 

This  research  was  conducted  as  a  four-phase  effort. 

•  Phase  1 :  State  of  the  Art  Review  and  Assessment:  This  phase  included  a  review  of 
related  publications,  research,  and  test  results.  Candidate  theories,  models,  test  methods  and 
historically  significant  field  test  data  were  identified  in  this  review.  Particular  emphasis  was 
placed  on  a  model  that  was  relatively  simple  to  calibrate  with  the  capability  to  capture  the 
critical  response  features  of  granular  material  behavior. 

•  Phase  2:  Model  Integration:  In  this  phase,  a  candidate  constitutive  model  was 
implemented  as  a  user  defined  material  model  in  the  ABAQUS  General  Purpose  Finite 
Element  Code. 

•  Phase  3:  Model  Calibration:  In  this  phase,  the  granular  material  response  model  was 
calibrated  with  laboratory  test  data.  The  testing  requirements  were  a  function  of  the  type  of 
model  selected  in  Phase  1 .  Historical  test  data  was  acquired  and  new  tests  were  conducted 
where  necessary  to  define  material  properties  for  unbound  granular  pavement  materials. 
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•  Phase  4:  Model  Verification,  Evaluation  and  Documentation:  In  this  phase,  the  newly 
calibrated  model  was  exercised  against  laboratory  test  data  and  selected  historical  field 
pavement  system  response  data  to  assess  the  predictive  suitability  of  the  model  (Webster, 

1 993).  The  ABAQUS  finite  element  code  was  used  to  make  these  predictions.  The  strengths 
and  weaknesses  of  the  response  model  and  calibration  parameter  relationships  were  evaluated 
and  documented. 
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CHAPTER  2:  PROBLEM  STATEMENT 


Classically,  the  flexible  pavement  used  in  military  airfields  consists  of  a  thin  asphalt 
concrete  (AC)  surface  to  provide  a  high-quality  waterproof  surface,  and  relatively  thick  layers 
of  granular  base  and  subbase  down  to  the  subgrade.  These  thick  granular  layers  are  used  to 
reduce  the  stresses  applied  by  aircraft  traffic  on  the  pavement  surface.  A  typical  pavement  of 
this  type  is  shown  in  Figure  2. 1 . 

The  magnitude  and  frequency  of  loading  in  airfield  pavements  are  very  different  from 
typical  highway  pavements.  The  magnitudes  of  aircraft  loadings  are  much  greater  than  the 
loads  seen  in  highways  as  shown  in  Figure  2.2.  The  amount  of  load  repetitions  applied  to 
airfield  pavement  is  several  orders  of  magnitude  less  than  that  seen  in  highways.  A  high- 
volume  highway  may  experience  60  million  equivalent  single  axle  loads  (ESAL),  while  a  high 
volume  airfield  may  only  experience  250,000  aircraft  coverages  in  a  20-year  period.  These 
differences  led  to  a  divergence  in  the  research  focus  between  the  airfield  and  highway 
pavement  communities.  The  major  focus  of  research  into  highway  flexible  pavement  design 
has  been  in  the  area  of  viscous  fatigue  modeling  of  asphalt  concrete.  The  airfield  pavement 
community  has  been  required  to  broaden  the  focus  of  analytical  research  to  include  the  AC 
and  all  supporting  layers  (Ahlvin,  1991). 

The  granular  base  and  subbase  have  always  posed  the  most  difficult  analytical 
problem  in  traditional  airfield  pavement  design  methodologies.  For  this  reason,  the  granular 
layers  have  never  been  treated  explicitly  in  design  as  have  the  AC  layer  and  subgrade  layer, 
which  have  used  predictive  models  for  cracking  in  the  AC  and  rutting  in  the  subgrade  as  a 
function  of  linear-elastic  strain  and  material  properties.  Instead  these  granular  layers  were 
carefully  specified  in  terms  of  gradation,  plasticity,  and  in-situ  density  to  minimize 
deformation  under  traffic.  However,  in  order  to  eventually  develop  theoretical  methods  to 
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predict  performance  of  the  pavement,  sound  methodologies  must  be  developed  that  will 
predict  plastic  deformation  within  these  granular  layers. 

The  structural  components  of  flexible  pavements  are  highly  nonlinear-elastic  plastic 
materials.  With  heavy  loading,  such  as  may  be  encountered  with  many  roads  and  airfields, 
the  nonlinear  response  of  pavement  materials  should  be  considered  when  predicting  pavement 
performance.  The  advances  made  in  computational  mechanics  have  created  new  tools,  such 
as  the  newer  generation  finite  element  codes,  for  this  type  of  problem.  The  beauty  of  the  finite 
element  method  is  that  it  can  incorporate  both  features  and  handle  arbitrary  geometries. 
Theoretically  rigorous  material  models  may  be  implemented  within  many  of  the  general- 
purpose  finite  element  computer  programs  available  today.  In  order  to  apply  these  material 
models,  mechanical  response  data  is  required  to  calibrate  the  necessary  model  parameters. 

The  essential  features  of  pavement  response  that  are  required  from  any  constitutive 
model  include  non-linear  elastic  response,  permanent  or  plastic  deformations  after  yield, 
cyclic  loading,  strain  softening/hardening,  and  shear  dilatancy.  This  research  addresses  the 
inadequacies  of  present  design  and  analysis  procedures  as  related  to  prediction  of  the  response 
of  granular  pavement  layers  subjected  to  aircraft  loads. 
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Figure  2. 1 .  Typical  flexible  pavement  configuration 
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Figure  2.2.  Comparison  of  aircraft  and  truck  (18K-ESAL)  loadings 
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CHAPTER  3:  LITERATURE  REVIEW 


MODEL  REQUIREMENTS 

Typical  rational  design  procedures  couple  theoretical  response  models  that  predict 
traffic  induced  stresses,  strains,  and  deflections  with  damage  models  for  fatigue  cracking  and 
pavement  rutting.  The  various  layers  in  a  pavement  system  are  characterized  by  their 
engineering  properties  and  the  structural  design  is  subsequently  based  upon  limiting  stresses, 
strains,  or  deflections  computed  at  certain  critical  locations  in  the  pavement  structure.  The 
procedures  use  an  iterative  process,  which  involves  theoretical  response  analysis,  material 
characterization,  distress  prediction,  and  adjustment  factors.  Several  rational  (mechanistic) 
pavement  design  procedures  have  been  introduced  into  design  over  the  past  years.  The 
development  of  the  procedures  is  summarized  in  the  Proceedings  of  the  International 
Conferences  on  the  Structural  Design  of  Asphalt  Pavements  (University  of  Michigan,  1962, 
1967, 1972,  1977,  1982,  1987,  and  International  Society  for  Asphalt  Pavements,  1992). 

ELASTICITY  MODELS 

Elasticity  models  can  be  divided  into  three  distinct  classes:  (1)  Cauchy  elasticity,  (2) 
hyperelasticity,  and  (3)  hypoelasticity.  Cauchy  elasticity  is  based  on  total  stress,  and  the 
current  stress  is  depends  only  on  the  current  strain.  Cauchy  elasticity  models  are  also 
reversible.  Hyperelasticity  is  a  total  stress  model  where  the  current  stress  depends  only  on  the 
current  strain.  In  addition,  hyperelasticity  models  are  based  on  the  principal  of  virtual  work  to 
insure  compliance  with  the  first  law  of  thermodynamics.  Hypoelastic  models  are  incremental 
stress  models  that  are  incrementally  reversible.  The  current  state  of  stress  ins  dependent  on 
the  stress  and  strain  path  followed.  Each  of  these  three  classes  of  models  are  addressed  in  the 
following  sections. 
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Many  theoretical  response  models  treat  a  pavement  system  as  a  layered,  linear  elastic 
system.  For  these  type  analyses,  load  associated  responses  are  governed  by  the  magnitude 
and  geometry  of  the  applied  loads,  and  the  thickness,  elastic  modulus,  and  Poisson’s  ratio  of 
the  individual  layers.  In  these  analyses,  each  layer  is  completely  characterized  by  the  elastic 
modulus  and  the  Poisson's  ratio.  Previous  research,  however,  has  shown  that  both  the  resilient 
(elastic)  and  permanent  deformation  behavior  of  granular  paving  materials  are  extremely 
complex,  depending  on  material  characteristics,  drainage,  and  loading  conditions.  The 
inability  of  layered  elastic  theory  to  account  for  stress  dependent  material  properties  has 
become  an  area  of  serious  concern. 

The  stress  dependency  of  the  resilient  modulus  and  Poisson’s  ratio  of  granular 
materials  have  been  examined  and  evaluated  by  a  number  of  researchers  (Seed,  et  al.,  1967, 
Hicks  and  Finn,  1970,  Hicks  and  Monismith,  1971,  and  Rada  and  Witczak,  1981).  The 
research  indicates  that  the  deformation  response  of  granular  materials  is  highly  stress 
dependent,  and  that  the  resilient  modulus  increases  with  increasing  confining  stresses.  The 
constitutive  relationship  developed  in  these  studies,  called  the  bulk  stress  model,  expresses  the 
resilient  modulus  as  a  function  of  the  bulk  stress,  using  Equation  3.1. 

Er=KOh  (3.1) 

where  Er  =  resilient  modulus 

2  =  bulk  stress  (Oi+2<t>3) 
kj,k2  =  regression  coefficients 

Also,  in  some  of  the  same  research,  Poisson’s  ratio  was  determined  to  increase  with 
increasing  values  of  the  principal  stress  ratio.  Historically,  the  bulk  stress  model,  above,  has 
been  combined  with  a  “constant”  Poisson's  ratio  and  used  widely  as  the  constitutive  model  of 
granular  materials  for  pavement  design. 
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In  continued  research  with  granular  materials.  May  and  Witczak  (1981)  and  Uzan 
(1985),  concluded  that  measured  and  predicted  responses  using  the  bulk  stress  model  did  not 
sufficiently  agree,  and  that  not  only  the  stress  state,  but  also  the  magnitude  of  the  induced 
shear  strains  influenced  the  resilient  modulus.  As  a  result,  Witczak  and  Uzan  (1988) 
improved  on  previous  relationships  and  introduced  the  octahedral  shear  stress  term  to  the 
determination  of  the  resilient  modulus.  In  addition  to  the  octahedral  shear  stress,  which  is  an 
invariant  shear  stress  term  for  three-dimensional  analysis,  they  also  made  the  equation 
dimensionally  correct  by  normalizing  the  bulk  and  shear  stress  terms  using  atmospheric 
pressure.  Equation  3.2  presents  the  modified  equation,  called  the  universal  model  because  it 
is  applicable  to  both  granular  and  cohesive  soils. 


xoct  =  octahedral  shear  stress 
PA  =  atmospheric  pressure 
ki,k2,k3  =  regression  coefficients 

The  use  of  Equations  3.2  and  3.3  above  have  been  recommended  as  the  constitutive 
model  for  unbound  pavement  material  layers  in  highway  performance  models  developed  at 
Texas  A&M  (Lytton,  et  al.,  1993).  The  five  material  parameters  in  these  models  are 
determined  using  nonlinear  regression  analysis  of  data  from  repeated  load  triaxial  tests. 

The  above  nonlinear  Cauchy  elastic  models  are  modifications,  or  a  simple  extension, 
of  the  generalized  form  of  Hooke’s  law,  and  use  secant  moduli  determined  from  the  stress  or 
strain  invariants;  thereby  accounting  for  confinement  effects.  These  Cauchy  elastic  models 
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are  total  stress  models  in  which  the  current  stress  depends  only  on  the  current  strain  and  the 
state  of  stress  is  path  independent.  An  essential  advantage  in  the  use  of  these  models  is  that 
the  model  parameters  have  physical  significance.  These  models  were  evaluated  by  Bonaquist 
(1996)  who  concluded  two  limitations.  First,  the  Cauchy  elastic  models  can  not  account  for 
the  volume  changes  which  result  from  the  application  of  shear  stresses,  because  they  are 
based  upon  Hooke’s  law  and,  therefore,  can  not  model  either  plastic  responses  or  dilation. 
Second,  nonlinear  elastic  models  may  violate  the  first  law  of  thermodynamics  and  generate 
energy  along  certain  cyclic  stress  paths  (Chen  and  Saleeb,  1982),  because  the  secant  moduli 
are  arbitrarily  selected. 

To  mitigate  the  problems  associated  with  violating  the  first  law  of  thermodynamics, 
hyperelastic  constitutive  relationships  have  been  developed  based  upon  the  principle  of 
conservation  of  energy  during  the  loading  and  unloading  of  an  elastic  body  (Lade  and  Nelson, 
1987,  Chen  and  Mizuno,  1990,  and  Uzan  et  al.,  1992).  Like  the  previously  mentioned  elastic 
models,  these  relationships  are  also  not  dependent  upon  the  stress  or  strain  history,  and  the 
stress-strain  behavior  is  both  reversible  and  path  independent.  Hyperelastic  models,  however, 
are  typically  higher  order  equations  with  a  large  number  of  regression  coefficients,  or  fitting 
parameters.  As  the  order  increases,  the  number  of  parameters  increases,  and  subsequently  the 
difficulty  in  performing  suitable  laboratory  tests  to  evaluate  the  parameters. 

One  of  the  more  straightforward  hyperelastic  models  is  that  proposed  by  Uzan  et  al. 
(1992),  as  a  part  of  the  Strategic  Highway  Research  Program.  They  assumed  both  the  resilient 
modulus  and  Poisson’s  ratio  to  be  stress  dependent  and  developed  the  following  stress 
dependent  relationship  for  Poisson’s  ratio  using  the  principle  of  conservation  of  energy  and 
the  universal  model  presented  in  Equation  3.2.  The  basic  form  of  this  non-linear  hyperelastic 
model  for  a  variable  Poisson’s  ration  is  given  in  Equation  3.3. 
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where  vs  =  secant  Poisson’s  ratio 

II  =  1st  stress  invariant  =  01+02+03 
J2  =  2nd  deviatoric  stress  invariant 
=  1/6  [(ai-02)2+(o2-a3)2+(o3-oi)2] 

01,02,03  =  principal  stress 

Bv(iJ)  =  Incomplete  Beta  function 
k2,k3,k4,k5  =  regression  coefficients 
Hyperelastic  models  are  total  stress  models,  which  satisfy  the  first  law  of 
thermodynamics,  account  for  nonlinearity,  confinement,  dilation,  and  can  be  used  to  model 
cyclic  loading.  Cyclic  loading  and  unloading,  however,  must  follow  the  same  path,  since  the 
current  stress  depends  on  the  current  strain.  The  primary  disadvantage  to  most  high  order 
hyperelastic  models  is  that  they  do  not  include  plastic  response  and  that  many  of  the  fitting 
parameters  have  no  physical  significance,  and  consequently  testing  to  evaluate  these 
parameters  is  frequently  complicated. 

A  third  type  elasticity  model  is  the  hypoelastic  constitutive  model,  which  addresses 
the  fact  that  in  many  materials,  including  granular  materials,  the  stress-strain  behavior  is  path 
dependent  and  the  response  is  not  necessarily  reversible.  The  hypoelastic  formulation  is  an 
incremental  constitutive  relationship,  with  the  behavior  determined  in  small  increments  of 
stress,  rather  than  for  the  entire  applied  stress.  The  current  state  of  stress  of  a  material 
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depends  upon  the  current  state  of  strain,  as  well  as,  the  stress  path  followed  to  reach  the 
current  state.  Like  the  hyperelastic  models,  the  true  hypoelastic  models  account  for 
nonlinearity,  confinement  effects,  and  dilation.  Unfortunately,  also  like  the  hyperelastic 
models,  many  of  the  true  hypoelastic  models  are  higher  order  formulations  which  result  in 
greater  numerical  complexity  and  a  large  number  of  material  fitting  parameters  which  have  no 
physical  significance,  or  interpretation. 

There  are,  however,  some  simpler  hypoelastic  models,  which  have  been  developed 
from  an  incremental  form  of  the  generalized  Hooke's  law.  For  these  models,  the  resilient 
modulus  is  replaced  with  variable,  or  incremental,  tangent  moduli,  which  are  functions  of  the 
stress  or  strain  invariants.  The  models  are  path  dependent  and  a  large  variety  of  nonlinear 
material  behavior  can  be  modeled.  While  these  models  lack  rigorous  theory  and  can  not 
include  dilation,  they  are  relatively  simple  and  the  model  parameters  do  have  physical 
significance.  Three  of  these  simplified  hypoelastic  models  are  presented  below.  Duncan  and 
Chang  (1970)  presented  the  following  model  shown  in  Equation  3.4. 


E  _T1  ^(l-sinOXtr,-^)!2 

'  _  2(ccosO  +  (T3sinO)  a 


(3.4  a) 


where: 


Et  =  tangent  Young’s  Modulus 
vt  =  tangent  Poisson’s  ratio 
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=  principal  stresses 
=  atmospheric  pressure 
=  cohesion 

=  angle  of  internal  friction 
=  material  constants 


Duncan  et  al.,  1978,  presented  the  following  as  in  Equation  3.5  where  the  elastic 
constants  are  functions  of  the  current  stress  state  and  the  Mohr  Coulomb  yield  surface 
location. 


E, 

Kt 


=  kbPa 


Rf  (1  -  sin  ®)(<7,  -  cr3 ) 
2 (c  cos®  +  <r3  sin®) 


(3.5  a) 


(3.5  b) 


where: 

E,  =  tangent  Young’s  Modulus 
IQ  =  tangent  bulk  modulus 
0^03  =  principal  stresses 
pa  =  atmospheric  pressure 
c  =  cohesion 

<p  =  angle  of  internal  friction 
k,n,kb,m  =  material  constants 
Rf  =  failure  ratio 

Domaschuk  and  Wade  (1969)  presented  the  following  relationship  shown  in  Equation 
3.6  ,  where  the  bulk  and  shear  modulus  constants  are  related  to  the  current  stress  state.  The 
octahedral  shear  and  normal  stresses  are  used  to  determine  the  elastic  constants. 
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K,  =K0  =mam 


(3.6  a) 


G,  =G0(l-bToc,)2  (3.6  b) 

where: 

K,  =  tangent  bulk  modulus 
G,  =  tangent  shear  modulus 
c  =  octahedral  normal  stress 
to,.,  =  octahedral  normal  stress 
Kq  =  initial  bulk  modulus 
G0  =  initial  shear  modulus 
b,m  =  material  constants 

In  the  use  of  hypoelastic  models,  initial  conditions  must  be  specified  since  the  stress- 
strain  behavior  of  the  materials  will  be  dictated  depending  upon  the  initial  starting  point.  With 
the  specification  of  loading  and  unloading  criteria,  these  models  can  be  used  to  model  the 
plastic  behavior  of  some  granular  materials.  These  models  have  seen  little  use  in  pavement 
analysis  except  for  limited  applications  in  nondestructive  pavement  testing  and  pavement 
thickness  design  to  resist  fatigue  cracking. 

PLASTICITY  MODELS 

Plasticity  models  characterize  the  plastic  deformation  behavior  of  soils  under  cyclic 
loading  and  are  particularly  useful  in  modeling  earthquake  responses.  Since  these  models 
predict  responses  to  cyclic  loading  their  benefits  in  performing  pavement  rutting  analyses  are 
obvious. 

The  first  type  models  considered  here  are  the  variable  modulus  models.  These 
models  are  based  upon  the  deformation  theory  of  plasticity,  are  relatively  simple  formulations 
derived  from  the  theory  of  elasticity.  The  incremental  nonlinear  elastic  constitutive  models 
presented  previously  in  are  frequently  used  with  the  variable  modulus  models  to  describe 
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permanent  deformation  under  cyclic  loading.  With  these  models,  different  tangent  moduli  are 
selected  (prescribed)  for  the  loading,  unloading,  and  reloading  conditions.  It  is  common 
practice  to  assume  that  the  unloading  and  reloading  moduli  are  equal  to  the  initial  tangent 
modulus  on  loading.  With  these  assumptions,  greater  deformations  occur  on  loading  than 
unloading,  and  as  a  result,  cyclic  loading  produces  permanent  deformations. 

Since  the  incremental  nonlinear  constitutive  models  presented  are  based  upon  the 
generalized  Hooke’s  law,  the  stress  invariants  for  octahedral  normal  stress  and  octahedral  shear 
stress  are  normally  used  to  define  volumetric  and  shear  loading  conditions.  If  more  complex 
constitutive  models  like  the  hyperelastic  or  hypoelastic  models  are  selected  for  use,  the 
distinction  between  loading  and  unloading  must  be  accomplished  with  the  use  of  an  energy 
density  function.  For  these  models,  loading  represents  positive  work,  while  unloading 
represents  negative  work. 

The  primary  advantage  to  variable  modulus  models  are  that  they  are  a  conceptually 
and  computationally  simple  formulation  and  a  logical  extension  to  the  elasticity  models  based 
upon  incremental  forms  of  Hooke’s  law,  presented  previously.  In  addition,  the  model 
parameters  used  in  the  models  have  physical  significance  and  interpretation.  A  disadvantage 
to  the  use  of  these  models  is  that  since  they  are  based  on  incremental  forms  of  Hooke’s  law, 
they  can  not  account  for  shear  dilation.  Another  disadvantage  is  that  these  models  violate 
continuity  for  the  neutral  loading  condition;  when  the  loading  function  is  equal  to  zero.  For 
this  condition,  either  loading  or  unloading  behavior  (and  moduli)  can  be  assumed. 

A  theoretically  rigorous  formulation  for  plasticity  has  been  developed  based  upon 
flow  theory.  Constitutive  models  based  upon  the  flow  theory  of  plasticity  are  incremental  and 
extend  the  elastic  stress-strain  relationships  into  the  plastic  range.  The  total  strain  is  the 
summation  of  the  reversible  elastic  strains  and  the  irreversible  plastic  strains.  Here  again,  an 
incremental  form  of  Hooke’s  law  using  elastic  moduli  (Young’s  modulus,  Lame’s  constants, 
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etc.),  that  are  functions  of  the  stress  or  strain  invariants  determines  the  elastic  strains.  The 
plastic  strains  are  functions  of  the  current  states  of  stress  and  strain,  and  the  incremental  stress 
gradient.  Yield  functions  are  introduced  in  flow  theory  to  differentiate  between  the  elastic  and 
plastic  states  (Chen  and  Mizuno,  1990)  (Salami,  1994). 

Yield  Functions 

A  yield  function  in  flow  theory  differentiates  between  elastic  and  plastic  behavior. 
Yield  functions  mathematically  describe  a  surface,  within  which  purely  elastic  recoverable 
deformations  or  strains  occur  and  along  which  purely  plastic  deformations  occur. 

Intersections  of  the  stress  path  with  the  yield  surface  result  in  both  elastic  and  plastic 
deformations.  Yield  functions  have  been  commonly  used  in  many  civil  and  geotechnical 
engineering  applications  to  describe  plastic  behavior  of  soils  and  other  construction  materials. 
Much  of  the  response  requirements  in  traditional  geotechnical  applications  require  monotonic 
loading  capabilities  only  (Chen  and  Mizuno,  1990),  while  pavements  applications  are  strongly 
tied  to  cyclic  response.  Five  of  the  more  commonly  applied  yield  functions  for  geotechnical 
materials  are  presented  below.  The  yield  functions  are  generally  expressed  in  terms  of  stress 
invariants  in  a  principal  stress  space  to  simplify  the  comparison  of  one  surface  to  another. 

The  basic  parameters  used  in  the  formulations  of  these  yield  functions  are  given  in  Table  3.1. 
Although  only  five  yield  functions  are  presented,  these  are  typical  the  large  number  of 
theories  that  have  been  proposed  over  the  last  40  years  of  geotechnical  engineering  history. 

Table  3.1.  Parameters  Used  in  Yield  Functions 
0  =  Lode  angle 

I,  =  first  invariant  of  the  stress  tensor 

J2  =  second  invariant  of  the  stress  deviator  tensor 

I3  =  third  invariant  of  the  stress  tensor 

J3  =  third  invariant  of  the  stress  deviator  tensor 

c  =  cohesion 

M=  angle  of  internal  friction 
k,a  =  material  constants 
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HYDROSTATIC 


Figure  3.1.  Coulomb  yield  function  and  surface 
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HYDROSTATIC 


4 J\  -  21 J I  -  36k1  J;  +  96k‘J,  -64k‘  =  0 

Figure  3.2.  Tresca  yield  function  and  surface 
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HYDROSTATIC 


alx  =  -  k  -  0 


Figure  3.3.  Drucker-Prager  yield  function  and  surface 


HYDROSTATIC 


J1-k 2  =  0 


Figure  3.4.  Von  Mises  yield  function  and  surface 
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The  Coulomb  yield  function  shown  in  Figure  3.1  is  a  three  dimensional  generalization 
of  the  well-known  Coulomb  failure  criterion  from  soil  mechanics.  This  yield  function  reduces 
to  the  Tresca  yield  function  for  the  case  of  frictionless  materials,  i.e.,  <f»  =  0.  While  both  of 
these  criteria  are  conceptually  simple,  they  both  have  singularities  at  the  comers  of  the 
hexagonal  shapes,  as  illustrated  in  Figures  3.1  and  3.2. 

These  singularities  are  avoided,  however,  with  use  of  the  Drucker-Prager  and  Lade- 
Duncan  yield  functions,  which  are  approximations  that  use  a  smooth  function  shown  in 
Figures  3.3  and  3.5.  In  addition,  the  Von  Mises  yield  function  (Figure  3.4)  is  a  smooth 
approximation  of  the  Tresca  yield  function  for  saturated  cohesive  soils  (frictionless  soils). 
While  both  the  Von  Mises  and  Drucker-Prager  smooth  yield  functions  neglect  the  effect  of  the 
third  stress  invariant,  the  Lade-Duncan  yield  function  includes  this  effect  in  the  approximation 
of  the  Coulomb  yield  function.  The  only  disadvantage  to  the  Lade-Duncan  approximation  is 
that  it  requires  multiaxial  testing  to  determine  the  material  coefficients;  otherwise,  it  is 
excellent  for  general  three-dimensional  analyses. 

Flow  Rules 

Once  a  yield  function  is  selected  to  differentiate  between  elastic  and  plastic  behavior 
of  granular  materials  analyzed  using  flow  theory,  a  flow  rule  is  needed  to  specify  the 
incremental  stress-strain  relationships  in  the  plastic  region.  The  flow  rule  specifies  the 
relationship  between  the  incremental  plastic  strains  and  the  current  state  of  stress  for  yielded 
materials  subjected  to  additional  loading.  The  flow  rule  states  that  the  direction  of  the  plastic 
strain  increment  is  normal  to  the  plastic  potential  function  at  the  current  state  of  stress 
(Bonaquist,  1996).  Equation  3.7  presents  the  general  mathematical  form  for  a  plasticity  flow 


where: 

de^j  =  plastic  strain  increment 
dA,  =  proportionality  constant  which  is  a  function  of  stress 
g  =  plastic  potential  function 
cjjj  =  current  state  of  stress 

If  the  plastic  potential  function  is  coincident  with  the  yield  function,  the  rule  is  said  to 
be  an  associated  flow  rule.  On  the  other  hand,  if  the  plastic  potential  function  is  different 
from  the  yield  function,  the  rule  is  called  a  non-associated  flow  rule. 

With  only  a  yield  function  and  flow  rule,  the  behavior  of  elastic  -  perfectly  plastic 
materials  can  be  modeled.  In  a  perfectly  plastic  material,  continued  loading  results  in  an 
increase  in  strain  with  no  increase  in  stress.  Soils  and  granular  materials,  however,  are  known 
to  exhibit  strain  hardening  or  strain  softening  with  continued  loading  in  the  plastic  region. 

Hardening  Rules 

Hardening  rules  have,  therefore,  been  developed  to  model  the  strain  hardening  and 
strain  softening  behaviors  of  soils  analyzed  using  flow  theory.  A  hardening  rule  permits,  and 
specifies,  a  movement  of  the  yield  function  in  stress  space  for  various  stress  increments.  An 
initial  yield  surface  (function)  is  specified,  and  once  the  stress  path  reaches  the  yield  surface 
subsequent  stress  increments  can,  and  normally  do,  result  in  the  generation  of  a  new  yield 
surface.  If  the  yield  surface  is  expanding,  hardening  behavior  is  said  to  be  occurring,  and  if 
the  yield  surface  contracts,  strain  softening  is  being  exhibited.  Stresses  within  the  yield 
surface  generate  elastic  responses,  while  stresses  that  intersect  the  yield  surface  result  in  a 
plastic  response. 

A  variety  of  hardening  rules  has  been  developed  to  model  the  behavior  of 
geotechnical  materials.  Most  plasticity  models  use  an  isotropic  hardening  rule,  which 
assumes  that  the  yield  surface  either  expands  or  contracts  uniformly  as  plastic  strains  occur. 
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If  a  yield  surface  translates  in  stress  space  as  a  rigid  body  and  retains  its  original  size  and 
shape,  a  kinematic  hardening  rule  is  applicable.  Complex  mixed  hardening  rules  are  possible 
which  allow  both  translation  and  expansion  or  contraction  of  the  yield  surface  as  plastic 
strains  occur.  One  of  the  better  known  applications  of  hardening  rules  is  the  cap  model, 
shown  schematically  in  Figure  1,  which  was  developed  specifically  for  geotechnical  materials. 

Cap  models  specify  a  failure  envelope,  above  which  plastic  behavior  occurs,  and  a 
strain-hardening  cap.  The  failure  envelope  is  typically  based  upon  one  of  the  yield  functions 
presented  earlier,  such  as  the  Drucker-Prager  or  the  Lade-Duncan  functions,  while  the  strain 
hardening  cap  can  be  modeled  with  a  variety  of  assumptions  -  an  ellipsoid,  a  sphere,  or  a 
straight  line.  With  this  type  model,  elastic  behavior  is  expected  when  the  stress  path  is  within 
the  “yield  surfaces,”  plastic  behavior  occurs  when  the  stress  path  intersects  the  failure 
envelope,  and  strain  hardening  occurs  when  the  stress  path  intersects  the  cap.  Bonaquist 
(1996)  considered  some  of  the  attributes  of  several  cap  models,  which  use  isotropic  hardening 
rules  as  shown  in  Table  3.2. 

The  Hierarchical  Single  Surface  (HiSS)  model  by  Desai,  et  al.  (1986)  was 
investigated  in  depth  by  Bonaquist.  The  HiSS  model  is  not  a  true  cap  model  since  it 
approximates  a  cap  model  using  a  single  continuous  function  to  provide  yield  surfaces  similar 
to  the  cap  models.  This  model  represents  a  significant  simplification  of  the  cap  models 
because  the  single  continuous  function  that  includes  both  the  yield  and  ultimate  failure 
surfaces  eliminates  singularities  and  the  numerical  difficulties  associated  with  the  two 
functions  used  in  conventional  cap  models.  This  model  has  since  been  extended  to  include 
strain  softening.  Bonaquist  concluded  that  the  HiSS  model  should  be  pursued  as  a 
constitutive  model  for  granular  pavement  materials.  The  HiSS  model  was  calibrated  for  a 

t 

base  course  material  and  verified  with  a  stand-alone  driver  code.  It  was  not  implemented  in  a 
finite  element  code  at  that  time. 


24 


Table  3.2.  Cap  Models  with  Isotropic  Hardening  Rules 


Model 

Yield 

Surface 

Flow  Rule 

Hardening  Rule 

Elastic 

Response 

Drucker  et  al., 
1957 

Drucker-Prager  with 
spherical  cap 

Associated 

Isotropic  hardening 
of  both  cone  and  cap 

Linear 

Di  Maggio  and 
Sandler,  1971 

Modified  Drucker- 
Prager  with  elliptical 
cap 

Associated 

Isotropic  hardening 
and  softening  of  cap 

Linear 

Sandler,  et  al., 
1976 

Modified  Drucker- 
Prager  with  elliptical 
cap 

Associated 

Isotropic  hardening 
of  cap 

General, 

incremental 

nonlinear 

Lade,  1975 

Modified  Lade- 
Duncan  with 
spherical  cap 

Two 

component: 

non- 

associated 

cone, 

associated  cap 

Isotropic  hardening 
and  softening  of  cone 
and  cap 

Duncan  and 
Chang  (1970) 
incremental 
nonlinear 

Baladi  and 
Rohani,  1979 

Drucker-Prager  with 
elliptical  cap 

Associated 

Isotropic  hardening 
of  cap 

General, 

incremental 

nonlinear 

(HiSS)  Desai 
et  al.,  1986 

Similar  to  Lade- 
Duncan  with  curved 
cap 

Associated  or 

non- 

associated 

Isotropic  or 
anisotropic  hardening 
of  cap  and  cone 

General, 

incremental 

nonlinear 

Isotropic  hardening  rules,  such  as  those  employed  in  the  plasticity  models  based  upon 
flow  theory  above  were  originally  developed  for  monotonic  loading  conditions.  Under  cyclic 
loading  conditions  these  models  are  inadequate.  While  they  extend  the  yield  surfaces  during 
elastic-plastic  behavior,  they  behave  elastically  during  unloading  and  reloading,  as  long  as  the 
stress  path  remains  within  the  yield  surface.  As  a  consequence,  cyclic  loading  at  the  same 
stress  state  result  in  no  additional  permanent  deformations.  To  account  for  the  hysteresis  and 
incremental  permanent  deformations  that  occur  during  cyclic  loading,  several  cyclic  load 
hardening  rules  have  been  developed  as  modifications  to  isotropic  hardening  models  (Desai  et 
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al.,  1986).  These  models  typically  consist  of  a  series  of  nested  yield  surfaces,  which  translate 
during  loading,  sequentially  intercepting  and  providing  different  yield  functions,  depending 
upon  the  stress  path  and  the  state  of  stress.  If  the  translating  yield  surfaces  are  allowed  to 
expand  and  contract,  depending  upon  the  state  of  stress,  complex  anisotropic  hardening 
models  can  be  generated  which  are  capable  of  modeling  a  wide  range  of  cyclic  behavior  and 
hysteresis  in  soils. 

Constitutive  models  based  upon  the  flow  theory  of  plasticity  provide  theoretically 
rigorous  solutions  and  numerical  stability  is  guaranteed  for  many  conditions.  These  models 
account  for  shear  dilation  and  the  model  parameters,  which  can  be  determined  from 
conventional  laboratory  triaxial  tests,  have  physical  significance.  The  primary  disadvantage 
to  flow  theory  models  is  that  the  numerical  analyses,  even  though  stable,  are  relatively 
complex  due  to  the  nature  of  the  yield  functions  in  stress  space  (Salami,  1994). 

On  the  other  end  of  the  spectrum  from  the  variable  modulus  models,  and  the  flow 
plasticity  models  are  the  theoretically  rigorous  formulations  for  plasticity  based  upon 
endochronic  theory.  Endochronic  theory  uses  incremental  constitutive  equations  and  extends 
the  elastic  stress  strain  relationships  into  the  plastic  range.  In  fact,  inelastic  behavior  is 
assumed  to  occur  from  the  onset  of  loading.  The  constitutive  relationships  divide  the  material 
responses  into  deviatoric  and  volumetric  components.  The  plastic  responses  are  subsequently 
characterized  by  scalar  variables  (intrinsic  time)  which  are  measures  of  the  rearrangement  of 
grain  configurations  during  plastic  deformation,  and  either  strain  hardening  or  strain  softening 
of  the  material.  For  rate  independent  materials,  the  scalar  variables  (sometimes  referred  to  as 
internal  variables)  are  functions  of  the  strain  history  and  related  to  the  length  of  the  plastic 
strain  path. 

Using  endochronic  theory  Valanis  (1971)  and  Valanis  and  Read  (1987)  developed 
constitutive  laws  for  the  inelastic  behavior  of  concrete  sand,  and  clay.  The  strain  hardening 
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and  strain  softening  functions  within  these  constitutive  models  are  determined  by  curve  fitting 
experimental  data  using  functional  forms,  which  represent  the  effects  of  structural  changes 
within  a  material.  Typically,  extensive  laboratory  testing  is  required  in  the  fitting  of  model 
parameters. 

A  wide  range  of  material  behavior,  including  cyclic  loading,  can  be  modeled  through 
the  appropriate  selection  of  elastic  constants  and  strain  hardening/softening  functions 
(Bonaquist,  1996).  These  strain  hardening/softening  functions  can  become  quite  complex 
when  a  large  range  of  material  behavior  is  modeled. 

The  complexity  of  formulation  and  extensive  laboratory  testing  required  for  plasticity 
models  has  been  a  traditional  source  of  reluctance  on  the  part  of  pavement  designers  to  use 
these  models  in  pavement  analysis  and  design  procedures.  Most  pavement  design  and 
construction  agencies  are  limited  to  traditional  geotechnical  and  materials  testing  capabilities 
(Ulidtz,  1998). 

RECENT  DEVELOPMENTS 

Recent  studies  have  explored  the  applicability  of  a  simpler  class  of  constitutive 
models  for  soils  based  on  micromechanics.  These  models  are  constructed  by  superposing  or 
integrating  the  response  of  smaller  units,  either  micromechanical  or  simply  mechanisms  of 
yielding  in  particular  stress  sub-spaces.  Often,  concepts  of  plasticity  are  stated  at  the  level  of 
the  postulated  micromechanism  in  order  to  characterize  its  kinetics.  The  numerical 
implementation  of  such  models  is  rather  delicate  (Peters,  1983)(Peters,  1997)  (Homer,  1997) 
(Prevost  and  Popescu,  1996). 

More  duly  micromechanically-  based  models  have  also  been  proposed.  In  these 
models,  soil  is  viewed  as  an  assemblage  of  particles,  and  the  unit  micro-mechanism  response 
is  defined  at  the  duly  micromechanical  level  of  contact  forces  with  rolling  and  sliding 
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kinematics  among  the  particles,  and  given  macroscopic  counter  parts  by  proper  definitions 
and  averaging  procedures  (Ulidtz,  1998)  (Prevost  and  Popescu,  1996). 

Recent  studies  at  WES  have  identified  a  relatively  simple  constitutive  model 
formulation  for  soils  that  is  a  non-linear  elastic-plastic  formulation  for  a  continuum  based  on 
response  laws  that  come  from  micromechanics.  The  model  recently  developed  at  WES 
(Peters,  1983,  1997,  1998)  has  been  used  successfully  in  vehicle  mobility  and  earthquake 
analysis  efforts  and  shows  great  promise  for  implementation  and  application  to  the  pavements 
problem 

The  elastic-plastic  model  produces  the  essential  features  of  soil  behavior  under 
complex  loading  histories  without  the  difficult  analytical  and  numerical  procedures  required 
for  calibration  and  implementation  of  existing  models  with  similar  capabilities.  The  central 
concept  is  a  multi-mechanical  model  that  produces  the  behavior  of  an  internal  variable  model; 
particularly  those  derived  from  endochronic  plasticity  theory.  As  for  an  endochronic  model, 
the  material  is  idealized  by  mechanisms  acting  in  parallel.  The  simplicity  comes  from  making 
each  mechanism  an  elastic-perfectly-plastic  element  that  approximates  the  response  of  an 
endochronic  element  (Valanis,  1971).  The  coupling  among  the  elements  is  mathematically 
simpler  than  for  the  endochronic  model,  a  feature  designed  to  simplify  both  calibration  and 
numerical  integration.  The  details  captured  best  by  the  model  are  initial  stiffness,  yield/failure 
stress,  shear-induced  volume  changes,  and  hysteresis  produced  by  cyclic  loading. 

SUMMARY  OF  LITERATURE  REVIEW 

This  review  of  potential  constitutive  models  for  granular  materials  in  pavements  has 
considered  a  number  of  different  model  formulations.  This  is  not  intended  to  be  a 
comprehensive  presentation  of  constitutive  models  for  soils  but  a  review  of  those  models  that 
have  received  attention  for  the  pavements  industry  to  date.  The  elastic  models  evaluated 
consisted  of  Cauchy  elastic,  hyperelastic,  and  hypoelastic  models.  Within  the  Cauchy  elastic 
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models,  the  bulk  stress  and  universal  models  were  considered  in  detail.  These  two  models  are 
nonlinear  extensions  of  the  generalized  form  of  Hooke’s  law,  which  use  the  secant  moduli, 
determined  from  the  stress  or  strain  invariants.  The  hyperelastic  models  considered  are  total 
stress  models,  which  satisfy  the  first  law  of  thermodynamics,  and  do  not  generate  energy 
along  certain  cyclic  stress  paths.  The  hypoelastic  models  likewise  satisfy  the  first  law  of 
thermodynamics,  but  address  the  fact  that  in  granular  materials  the  stress-strain  behavior  is 
path  dependent,  and  the  response  is  not  necessarily  reversible. 

The  plastic  models  considered  in  this  review  addressed  formulations  based  upon 
plasticity  theory,  endochronic  theory,  micromechanical  theory,  and  the  WES  Multimechanical 
elastic-plastic  model.  A  flow  diagram  summarizing  the  evolution  of  these  plastic  models  is 
shown  in  Figure  3.  6.  Models  based  upon  the  deformation  theory  of  plasticity  represent 
extensions  of  incremental,  nonlinear  elastic  models,  and  extend  such  models  to  cover  both 
loading  and  unloading  behavior.  The  models  based  upon  endochronic  theory  use  no  loading 
criteria,  or  yield  surfaces,  and  elastic-plastic  response  is  assumed  from  the  beginning  of 
loading.  With  these  models,  a  scalar  internal  variable  called  intrinsic  time  is  used  to  account 
for  loading  history  and  the  stress  path.  Finally,  models  based  upon  flow  theory  of  plasticity 
were  considered.  These  models  extend  the  elastic  stress-strain  relationships  into  the  plastic 
region  with  the  use  of  a  yield  function  which  differentiates  between  elastic  and  elastic-plastic 
material  behavior.  The  yield  function  defines  a  surface  in  stress  space,  inside  of  which  elastic 
behavior  occurs,  and  on  and  outside  of  which  plastic  responses  can  be  expected.  A  flow  rule 
is  used  to  specify  the  incremental  stress-strain  relationships  that  occur  in  the  plastic  region, 
i.e.,  outside  the  yield  surface.  Strain  hardening  and  strain  softening  behaviors  are  modeled  by 
specifying  hardening  rules,  which  permit  movement  of  the  yield  surface  in  stress  space  for 
various  stress  increments. 
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The  use  of  elastic  models  for  pavement  systems  has  been  generally  restricted  to 
nondestructive  pavement  testing  and  pavement  thickness  design  to  resist  fatigue  cracking  and 
subgrade  rutting,  due  to  their  inability  to  adequately  model  cyclic  loading.  On  the  other  hand, 
plasticity  models,  with  their  ability  to  model  cyclic  loading  and  plastic  deformations,  are 
obviously  beneficial  in  modeling  rutting  behavior  and  permanent  deformations  in  pavement 
systems.  However,  each  plasticity  theory  has  certain  advantages  and  disadvantages  when  it 
comes  to  implementation  in  pavement  systems  modeling. 

Models  based  upon  deformation  theory  of  plasticity  are  direct  extensions  of  the 
incremental  forms  of  Hooke’s  law,  and  as  a  result  are  conceptually  straightforward  and 
computationally  simple.  In  addition,  the  model  parameters  used  in  deformation  models  have 
physical  and  engineering  significance.  On  the  other  hand,  these  models  can  not  account  for 
shear  dilation  and  violate  continuity  conditions  for  neutral  loading  conditions. 

Endochronic  models  use  relatively  straightforward  constitutive  relationships  and  use  a 
scalar  internal  variable  to  govern  inelastic  responses  and  account  for  strain  history.  Unlike 
deformation  theory  models,  endochronic  models  can  model  and  account  for  shear  dilation. 

The  theory  is  relatively  new,  however,  and  currently  only  limited  applications  have  been 
developed.  When  used,  model  parameters  have  physical/engineering  significance,  but  fitting 
of  model  parameters  requires  extensive  laboratory  materials  testing. 

Models  based  upon  flow  theory  of  plasticity  provide  theoretically  rigorous  solutions, 
and  numerical  stability  is  guaranteed  for  certain  conditions.  These  models  can,  also,  account 
for  shear  dilation  and  their  parameters  have  physical  and  engineering  significance.  Material 
behavior  is  divided  into  elastic  and  elastic-plastic  responses  by  yield  functions,  which  can  be 
relatively  complex  shapes  in  stress  space.  Incremental,  nonlinear  elastic  models  are  used 
inside  the  yield  surfaces  to  define  material  behavior,  while  flow  rules  and  hardening  rules  are 
use  to  define  the  response  on  and  outside  the  yield  surface. 
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To  date  plasticity  models  have  not  been  used  extensively  in  pavement  applications. 
This  fact  is  a  result  of  several  factors.  First,  their  primary  application  would  be  in  modeling 
rutting  and  permanent  deformations  in  pavement  systems,  which  typically  result  from  cyclic 
or  repeated  load  applications.  The  modeling  of  cyclic  or  repeated  load  applications  using 
plasticity  models  is  computationally  intensive,  requiring  the  dedication  of  significant 
computing  resources.  Next,  typical  values  of  model  parameters  for  most  common  paving 
materials  have  not  been  established,  and  can  not  be  derived  from  traditional  empirical 
characterization  tests  used  for  soil  and  aggregate  bases.  In  addition,  the  majority  of 
geotechnical  tests  performed  on  soils  and  aggregates  do  not  evaluate  the  effects  of  cyclic 
hardening  or  softening  of  the  materials. 

A  constitutive  model  that  can  capture  the  essential  behavior  of  pavement  materials 
under  service  environments  has  many  requirements  including  simplicity  of  calibration  and 
operation,  physical  significance  of  the  model  parameters,  and  the  ability  to  be  readily 
incorporated  into  analysis  codes.  The  WES  Multimechanical  model  possesses  all  of  these 
features  and  is  yet  untested  in  the  pavement  community,  and  its  application  to  pavement 
system  analysis  will  be  the  primary  focus  of  this  dissertation. 
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Development  Of  Plasticity 
Theory  Based  On  A  Yield 
Surface 

Isotropic  Hardening 
Kinematic  Hardening 
1950-60s 

Critical  State  Theory 
(Soil  Mechanics) 

1960s 

Multi-Surface  Models 
1970s 

Bounding  Surface  Models 
1980s 

Hypoplasticity 

1990s 


Irreversible 
Thermodynamics  and 
Internal  Variable  Theory 
1950s 


General  Treatment  of 
Viscoelasticity 
1960s 


Micro-Mechanical 

Application  of  Internal 

Theories 

Variable  Theory  to  Rate- 

Independent  Materials 

(Plasticity) 

1970s-Present 

(e.g.  Endochronic  Theory) 

1970s-Present 

Multi-Mechanism  Model 
(An  interpretation  of  internal 
variable  theory) 


Figure  3.6.  Flow  diagram  summarizing  the  evolution  of  plastic  constitutive  models  for  soils 
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CHAPTER  4:  SELECTION  AND  IMPLEMENTATION 

OF  MODEL 

CANDIDATE  CONSTITUTIVE  MODELS 

The  essential  features  of  pavement  response  that  are  required  from  any  constitutive 
model  include  non-linear  elastic  response,  permanent  or  plastic  deformation  after  yield,  cyclic 
loading,  strain  softening/hardening,  and  shear  dilatancy.  A  pavement  model  should  be  simple 
in  operation,  implementation  and  calibration.  The  model  must  be  executable  within  a  proven 
general  purpose  Finite  Element  code.  Of  the  general  classes  of  constitutive  theories  studied 
(Linear  Elastic,  Non-Linear  Elastic,  and  Plasticity),  only  those  theories  based  on  plasticity 
have  the  necessary  features  to  perform  adequately  as  a  model  for  granular  pavement  materials. 
A  summary  of  selected  models  discussed  in  Chapter  3  and  their  features  is  shown  in  Table  4.1. 

The  HiSS  model  by  Desai  was  thoroughly  investigated  by  Bonaquist  in  1996  at  the 
University  of  Maryland.  Although,  Bonaquist  concluded  that  it  shows  promise  as  a  potential 
model  for  granular  pavement  material,  the  HiSS  model  does  not  appear  to  have  the  simplicity 
of  calibration  and  implementation  desired  for  a  pavement  material  model. 

Many  engineers  in  the  pavement  industry  tasked  with  advanced  analysis  of  pavement 
behavior  will  use  a  commercial  general  purpose  FEM  like  ABAQUS  as  their  typical  analysis 
program  since  special  purpose  non-linear  FEM  programs  for  pavements  are  not  readily 
available.  The  Modified  Drucker-Prager  (DP)  is  recommended  by  ABAQUS  as  the  model  for 
use  in  modeling  granular  material  behavior.  The  DP  model  has  been  around  in  various  forms 
for  many  years  and  was  originally  developed  for  soils  with  much  lower  strength  than  the  base 
course  materials  under  investigation  here.  It  has  the  capability  to  capture  ultimate  failure/yield 
stress  for  a  wide  range  of  materials,  however  it  does  not  have  the  sophistication  required  to 
adequately  represent  the  complex  multi  -stage  yielding  seen  in  highly-compacted  granular 
materials.  Its  usefulness  for  this  effort  is  to  demonstrate  the  inadequacies  of  classical 
constitutive  models  that  one  would  find  in  an  FEM  code  like  ABAQUS. 
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The  WES  Multimechanical  elastic-plastic  model  produces  essential  features  of  soil 
behavior  without  the  difficult  analytical  and  numerical  procedures  required  for  calibration  and 
implementation  of  existing  models  with  similar  capabilities.  The  details  captured  best  by  the 
model  are  initial  stiffness,  yield/failure  stress,  shear-induced  volume  changes,  and  cyclic 
behavior.  The  WES  Multimechanical  model,  which  shows  high  potential  in  the  area  of 
granular  pavement  material  modeling,  its  calibration  requirements,  and  its  application  for 
constitutive  modeling  of  granular  pavement  materials  will  be  the  primary  focus  of  this 
research.  A  discussion  of  the  ABAQUS  Drucker-Prager  Cap  model  and  the  WES 
Multimechanical  model  follows. 


Table  4. 1 .  Critical  Features  of  Selected  Models  for  Unbound  Pavement  Materials 


Model 

Critical  Response  Features 

Linear 

Elastic 

Non-Linear 

Elastic 

Plastic 

Shear 

Dilation 

Cyclic 

Loading 

Selected  Reference 

Linear 

Elasticity 

X 

Barker  and  Gonzalez, 

1991 

Bulk  Stress 

(KThetaK) 

X 

X 

Hicks  and  Monismith, 

1971  Rada  and  Witczak, 

1981 

HiSS 

X 

X 

X 

X 

X 

Desai,  1986,  Bonaquist, 

1996 

WES 

X 

X 

X 

X 

X 

Homer,  1997,  Peters, 

1998,  Meade,  1997, 

1998 

Drucker- 

Prager 

X 

i 

X 

Limited 

Baladi  and  Rohani, 

1979,  ABAQUS  Theory 

Manual,  1998 
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DESCRIPTION  OF  ABAQUS 

ABAQUS  is  a  general-purpose  finite  element  program  developed  and  marketed  by 
Hibbitt,  Karlsson,  and  Sorensen,  Inc.  of  Pawtucket,  Rhode  Island.  ABAQUS  is  written  in 
transportable  FORTRAN,  although  the  input/output  routines  are  optimized  for  specific 
computer  systems.  The  source  code  for  ABAQUS,  not  available  to  the  user,  contains  about 
300,000  executable  statements. 

One  of  the  most  important  features  of  ABAQUS  is  its  use  of  the  library  concept  to 
create  different  models  by  combining  different  solution  procedures,  element  types,  and 
material  models.  The  analysis  module  consists  of  an  element  library,  a  material  library,  a 
procedure  library,  and  a  loading  library.  Selections  from  each  of  these  libraries  can  be  mixed 
and  matched  in  any  reasonable  way  to  create  a  finite  element  model. 

The  material  library  includes  linear  and  nonlinear  elasticity  models  as  well  as 
plasticity  and  viscoplasticity  formulations.  The  analysis  procedure  library  includes  static 
stress  analysis,  steady  state  and  transient  dynamic  analysis,  and  a  number  of  other  specialized 
procedures.  In  all  of  these  analysis  types,  time  is  used  as  the  index  for  incremental  solution 
techniques.  Time  is  a  purely  arbitrary  index  in  the  static  procedures  used  in  this  study. 

ABAQUS  DRUCKER-PRAGER  CAP  MODEL 

The  modified  Drucker-Prager  /  Cap  plasticity  is  intended  to  model  cohesive 
geological  materials  that  exhibit  pressure-dependent  yield,  such  as  soils  and  rocks.  It  is  based 
on  the  addition  of  a  cap  yield  surface  to  the  Drucker-Prager  plasticity,  which  provides  an 
inelastic  hardening  mechanism  to  account  for  plastic  deformation  and  helps  to  control  volume 
dilatancy  under  yielding.  The  ABAQUS  DP  model  provides  a  reasonable  response  to  large 
stress  reversals  in  the  cap  region  through  an  isotropic  hardening  rule;  however,  in  the  failure 
surface  region  the  response  is  reasonable  only  for  essentially  monotonic  loading  (ABAQUS, 
1998). 
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Yield  Surface 


The  addition  of  the  cap  yield  surface  to  the  Drucker-Prager  model  serves  two  main 
purposes:  it  bounds  the  yield  surface  in  hydrostatic  compression,  thus  providing  an  inelastic 
hardening  mechanism  to  represent  plastic  compaction.  The  addition  of  the  cap  also  helps  to 
control  volume  dilatancy  when  the  material  yields  in  shear  by  providing  softening  as  a 
function  of  the  inelastic  volume  increase  created  as  the  material  yields  on  the  Drucker-Prager 
shear  failure  surface. 

The  yield  surface  has  two  principal  segments:  a  pressure-dependent  Drucker-Prager 
shear  failure  segment  and  a  compression  cap  segment,  as  shown  in  Figure  4.1.  The  Drucker- 
Prager  failure  segment  is  a  perfectly  plastic  yield  surface  (no  hardening).  Plastic  flow  on  this 
segment  produces  inelastic  volume  increase  (dilation)  that  causes  the  cap  to  soften.  On  the 
cap  surface,  plastic  flow  causes  the  material  to  compact. 

Transition 


Figure  4.1.  ABAQUS  Drucker-Prager  model  yield  surface 
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Failure  Surface 


The  ABAQUS  Drucker-Prager  failure  surface  is  written  in  a  q  (principal  stress 
difference,  q=crrcr3)  versus p  (mean  normal  stress , p=(<Ji+2cr3)/3)  space  as: 

Fs  =q-ptanj3-d  =  0  (4.1) 

Where,  cr,  is  the  maximum  principal  stress,  cr3  is  the  minimum  principal  stress,  /?  represents 
the  angle  of  friction  in  the  q-p  plane,  and  d  is  the  cohesion. 

Cap  Yield  Surface 

The  cap  yield  surface  has  an  elliptical  shape  with  constant  eccentricity  in  q-p  plane 
and  also  includes  dependence  on  the  third  stress  invariant  in  the  deviatoric  plane.  The  cap 
surface  hardens  or  softens  as  a  function  of  the  volumetric  inelastic  strain.  The  ABAQUS 
Drucker-Prager  cap  yield  surface  Fc  and  transition  surface  Ft  is  written  as  : 


Fc=]\p~Pal  + 


Rq 


-|2 


|_(l  +  a-a/cos/?)J 


-R(d  +  pa  tan  /?)  =  0 


(4.2  a) 
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“  2 

Ip-pJ  + 
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1 

/I 

(d  +  pa  tan/?) 

|( 

_ 

^  COS  j3  j 

a(d  +  pa  tan  /?)  =  0 


(4.2  b) 


Where  R  is  a  parameter  that  controls  the  shape  of  the  cap,  a  is  a  cap  transition  factor,  and  pa 
is  an  evolution  parameter  that  represents  the  volumetric  inelastic  strain  driven 
hardening/softening.  The  pa  parameter  is  a  function  of  the  plastic  volumetric  strain  and 
volumetric  yield  stress  pb. 


Defining  Yield  and  Hardening  Parameters 

The  variables  d,  [},  R,  and  a  are  provided  by  the  user  to  define  the  shape  of  the  yield 
surface.  The  hardening  curve  specified  for  this  model  interprets  yielding  in  the  hydrostatic 
pressure  sense:  the  hydrostatic  pressure  yield  stress  is  defined  as  a  tabular  function  of  the 
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volumetric  inelastic  strain,  and,  if  desired,  a  function  of  temperature  and  other  predefined  field 
variables.  The  range  of  values  for  which  pb  is  defined  should  be  sufficient  to  include  all 
values  of  effective  pressure  stress  that  the  material  will  be  subjected  to  during  the  analysis. 

Plastic  Flow 

Plastic  flow  is  defined  by  a  flow  potential  that  is  associated  in  the  deviatoric  plane, 
associated  in  the  cap  region  in  the  meridional  plane,  and  nonassociated  in  the  failure  surface 
and  transition  regions  in  the  meridional  plane.  The  flow  potential  surface  is  made  up  of  an 
elliptical  portion  in  the  cap  region  that  is  identical  to  the  cap  yield  surface,  and  another 
elliptical  portion  in  the  failure  and  transition  regions  that  provides  the  nonassociated  flow 
component  in  the  model.  The  two  elliptical  portions  form  a  continuous  and  smooth  potential 
surface  (ABAQUS,  1998). 

Calibration 

At  least  three  experiments  are  required  to  calibrate  the  simplest  version  of  the  DP 
model:  a  hydrostatic  compression  test  and  two  triaxial  compression  tests  (more  than  two  tests 
are  useful  for  a  more  accurate  calibration).  A  more  detailed  discussion  of  the  tests  and 
procedures  used  for  calibration  is  given  in  Chapter  5. 

WES  MULTIMECHANICAL  CONSTITUTIVE  MODEL 

Background 

The  elastic-plastic  model  produces  the  essential  features  of  soil  behavior  under 
complex  loading  histories  without  the  difficult  analytical  and  numerical  procedures  required 
for  calibration  and  implementation  of  existing  models  with  similar  capabilities.  The  central 
concept  is  a  multi-mechanical  model  that  mimics  the  behavior  of  internal  variable  model, 
particularly  those  derived  from  endochronic  plasticity  theory.  As  for  an  endochronic  model, 
the  material  is  idealized  by  mechanisms  acting  in  parallel.  The  WES  model  uses  four 
mechanisms  in  its  current  form.  The  simplicity  comes  from  making  each  mechanism  an 
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elastic-perfectly-plastic  element  that  approximates  the  response  of  an  endochronic  element. 
The  coupling  among  the  elements  is  mathematically  simpler  than  for  the  endochronic  model,  a 
feature  designed  to  simplify  both  calibration  and  numerical  integration.  The  details  captured 
best  by  the  model  are  initial  stiffness,  yield/failure  stress,  shear-induced  volume  changes,  and 
hysteresis  produced  by  cyclic  loading. 

In  order  to  accomplish  the  objectives  of  this  research  the  WES  model  was 
implemented  in  two  distinct  forms.  A  PC-Compatible  stand-alone  version  and  an  ABAQUS 
User  Defined  Material  Model  Subroutine,  (UMAT).  The  stand-alone  model,  MVTEWER,  was 
used  to  provide  quick  feedback  during  the  iterative  calibration  process  for  the  WES  model.  A 
discussion  of  the  MVIEWER  program  is  presented  in  Appendix  G.  The  MVIEWER  was 
compiled  using  a  commercial  PC  compatible  FORTRAN  77  compiler.  Since  this  model  had 
originally  been  developed  for  use  on  a  PC  it  was  relatively  simple  to  take  the  model 
subroutines  and  add  a  constitutive  driver  program  to  produce  outputs  of  stress  and  strain  for  a 
given  stress  or  strain  path.  The  UMAT  was  programmed  in  FORTRAN  77  according  to  the 
guidelines  given  by  ABAQUS  for  development  and  implementation  of  a  user-defined  material 
model. 

General  Description 

The  elastic-plastic-perfectly-plastic  elements  act  in  parallel  by  making  the  total  strain 
common  to  all  mechanisms  as  represented  in  Figure  4.2.  Thus,  each  element  is 
computationally  independent  and  can  be  integrated  using  an  efficient  radial  return  procedure. 
The  total  stress  is  the  sum  of  the  component  stresses.  The  shear  and  hydrostatic  mechanisms 
are  independent  because  they  represent  different  deformation  mechanisms.  A  coupling  exists 
between  shear  and  hydrostatic  mechanisms  in  the  form  of  a  shear-dilatancy  law.  The  coupling 
law  imparts  a  plastic  hydrostatic  strain  increment  to  the  total  volumetric  strain  that  in 
proportion  to  the  total  plastic  shear  strain  produced  by  the  shear  mechanisms.  The  volumetric 
proportionality  constant  depends  on  the  shear  stress  to  hydrostatic  stress  ratio  in  a  manner 
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reminiscent  of  classical  Critical  State  Soil  Mechanics  (CSSM)  (Schofield,  C.  P.,  and  Wroth, 
D.  M.,  1968).  In  contrast  to  the  CSSM  unidirectional  dilatancy  law,  the  present  model  senses 
the  direction  of  shear  loading  and  correctly  predicts  the  magnitude  and  sign  of  plastic 
volumetric  strain  during  unloading. 

The  stresses  within  the  mechanism,  and  the  void  ratio  of  the  soil  describe  the  material 
state.  The  plastic  strains  are  thermodynamic  “forces”  that  retain  the  effects  of  the  stress 
history  of  the  material.  The  model  uses  three  groups  of  parameters:  stiffness  parameters, 
strength  parameters,  and  a  shear-volume  coupling  parameters.  (Meade,  1998)  (Peters,  1998) 

Stiffness  Parameters 

The  stiffness  parameters  are  shear  modulus  for  each  shear  mechanism  and  bulk 
modulus  for  each  hydrostatic  mechanism.  The  sum  of  the  stiffness  moduli  defines  the  initial 
elastic  stiffness  of  the  material.  By  distributing  the  moduli  among  the  mechanisms  according 
to  the  mechanism’s  yield  strength,  the  shape  of  the  stress-strain  curve  can  be  modeled. 


Figure  4.2.  Idealized  representation  of  WES  Multimechanical  model 
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Yield  Parameters 


The  strength  parameters  define  the  yield  stress  for  each  mechanism.  Each  mechanism 
acts  as  an  elastic-plastic  component  whereby  the  response  is  elastic  for  all  stress  increments 
within  the  surface  and  plastic  when  the  stress  point  lies  on  the  surface.  Stress  increments  that 
fall  outside  of  the  surface  are  scaled  back  to  the  surface. 

A  friction  parameter  and  cohesion  determine  the  limiting  shear  stress.  The  friction  is 
introduced  through  a  yield  law  of  the  form: 

non -r  M  (4.3) 

where  Qr  is  the  total  stress  for  mechanism  r  defined  as: 

Qr  =  Qrs+a '  (a+a)  (4.4) 

The  shear  component  Qr  s  is  determined  from  the  constitutive  response  of  the  mechanism. 

The  hydrostatic  component,  (cr+a)  is  distributed  from  the  total  hydrostatic  stress  and  cohesion 
in  proportion  to  the  distribution  factor  ar.  Thus  the  shear  mechanism  sees  the  hydrostatic 
stress  as  a  parameter.  The  function /is  chosen  to  represent  a  Mohr-Coulomb-like  yield 
surface  with  Y  being  the  limit  parameter  for  the  mechanism  that  is  scaled  to  the  friction  angle, 
<j>,  of  the  material. 

Yield  of  the  hydrostatic  mechanisms  is  scaled  by  a  reference  stress  that  depends  on 
void  ratio  by  the  law: 

Qk  =  VPe  (e)  (4.5) 

The  scale  factor  K  determines  the  limit  stress  of  hydrostatic  mechanism,  r.  The  reference 
stress,  Pe(e),  lies  on  the  virgin  loading  curve  at  the  point  corresponding  to  the  prevailing  void 
ratio,  e.  The  effect  of  void  ratio  on  shear  response  comes  through  the  dependence  of  shear 
yield  stress  on  the  hydrostatic  stress. 

Materials  possessing  cohesion  can  withstand  some  tensile  stresses.  The  tensile 
strength  is  accounted  for  by  applying  a  reduction  to  the  mean  stress  that  is  proportional  to  the 
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material  cohesion.  Each  mechanism  is  allocated  a  portion  of  the  tensile  strength  in  proportion 
to  the  amount  of  volumetric  stiffness  that  the  mechanism  contributes  to  the  overall  bulk 
stiffness  of  the  material. 

Shear-Volume  Coupling 

The  magnitude  of  the  shear  volume  coupling  is  controlled  by  two  parameters,  the 
ratio,  Mc,  of  shear  to  hydrostatic  stress  at  which  a  specimen  begins  to  dilate  in  a  monotonic 
loading  test  and  a  parameter,  y,  that  scales  the  dilatancy  rate  as  the  stress  ratio  becomes 
greater.  In  the  CCSM,  only  Mc  is  used  because  it  is  assumed  by  critical  state  theory  y=l .  The 
hydrostatic  strain  “seen”  by  the  hydrostatic  mechanisms  is  distinct  from  that  caused  by 
coupling  with  the  shear. 

Details  of  Calculations 

The  model  computation  is  strain  driven.  Given  the  current  internal  state  and  strain 
increment,  the  model  produces  an  updated  stress  state.  The  integration  procedure  is  explicit. 
First  the  response  for  each  shear  mechanism  is  computed.  This  computation  consists  of  (1) 
computing  an  elastic  “trial”  stress,  (2)  comparing  the  resulting  elastic  stress  to  the  yield  stress, 
and  (3)  if  beyond  yield,  scaling  back  along  a  radial  path  to  the  yield  surface.  The  elastic  strain 
associated  with  the  stress  increment  inside  the  yield  surface  is  subtracted  from  the  total  strain 
increment.  The  plastic  shear  strain  is  computed  as  the  difference  between  the  elastic  stress  and 
total  stress  divided  by  the  shear  modulus. 

Once  the  shear  response  is  computed,  the  plastic  shear  strain  is  used  to  compute  the 
volumetric  strain  that  results  from  shear-volume  coupling.  The  shear  strain  is  the  weighted 
sum  of  the  shear  strain  for  the  individual  mechanisms.  The  weighting  factor  for  a  mechanism 
is  the  ratio  of  its  shear  modulus  to  the  total  shear  modulus.  These  factors  add  up  to  one.  The 
dilatancy  strain  is  removed  from  the  total  volumetric  strain  to  produce  a  net  hydrostatic  strain 
that  is  used  for  the  computation  of  the  hydrostatic  response. 
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The  hydrostatic  response  is  similar  to  that  of  the  shear  stress.  A  trial  elastic  stress  is 
computed  which  is  compared  to  the  limit  hydrostatic  stress.  If  the  stress  exceeds  the  limit 
stress,  it  is  scaled  back  to  the  limit  value. 

Finally,  the  shear  stress  is  adjusted  to  account  for  a  reduction  in  hydrostatic  stress  due 
to  combined  effects  of  dilatancy  and  hydrostatic  strain.  The  adjustment  is  accomplished 
simply  by  setting  the  shear  strain  increment  set  to  zero  and  using  the  shear  computation 
described  previously.  Note  that  the  computation  for  the  shear  strain  treats  the  hydrostatic 
stress  as  a  parameter.  Thus  this  final  step  can  be  viewed  as  an  adjustment  to  account  for  a 
change  in  a  state  dependant  parameter. 

The  numerical  procedure  is  efficient,  without  iteration,  and  is  accurate.  It  does  not 
become  unstable  near  failure  and  its  efficiency  is  virtually  the  same  for  both  elastic  and  plastic 
conditions. 

Coding  Details 

The  model  has  been  implemented  in  the  finite  element  program  ABAQUS.  The 
ABAQUS  program  permits  the  user  to  write  a  subroutine  that  contains  a  user-defined 
constitutive  model  or  UMAT.  The  UMAT  was  written  in  FORTRAN  77  and  consists  of  one 
main  subroutine,  five  sub-task  subroutines,  and  two  functions.  A  separate  subroutine,  titled 
SDVINI,  was  written  to  initialize  solution-dependent  state  variables,  which  include  the  full 
stress  tensor  and  a  void  ratio.  The  FORTRAN  source  code  is  shown  in  Appendix  A. 

ABAQUS  Features 

ABAQUS  is  a  general-purpose  finite  element  program  licensed  from  Hibbitt,  Karlsson 
&  Sorenson,  Inc.  Version  5.8  of  ABAQUS  was  used.  The  program  permits  the  user  to 
employ  a  constitutive  model  of  one’s  choosing.  The  model  calculations  are  contained  in  the 
UMAT.  The  UMAT  author  must  conform  to  certain  conventions  to  enable  the  UMAT  to 
interface  properly  with  the  finite  element  solver.  The  user  may  specify  material  properties  that 
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are  entered  on  a  command  line  in  the  input  data  set.  The  user  may  use  solution-dependent 
state  variables  that  may  be  updated  within  the  UMAT. 

The  main  program  calling  the  UMAT  provides  stresses  at  the  start  of  a  loading  step, 
total  strains,  and  strain  increments  for  the  current  step.  The  UMAT  must  determine  the  stress 
increment  caused  by  the  strain  increment  and  update  both  the  stress  and  solution-dependent 
state  variables  at  end  of  the  step.  In  addition,  the  UMAT  must  provide  a  stress  gradient 
matrix,  a  Jacobian.  The  Jacobian  is  an  estimate  of  the  stiffness  at  the  current  material  state, 
which  uses  the  most  recent  stresses  that  were  in  equilibrium.  A  direct  strain  increment  is 
applied  to  each  direction  X  and  Y,  and  a  shear  strain  increment  is  applied  to  the  X-Y  plane. 
These  strains  are  applied  independently  and  the  stress  increment  produced  by  each  strain  is 
calculated.  The  ratio  of  stress  increment  to  the  strain  increment  is  used  as  an  estimate  of  the 
Jacobian. 

Material  Properties 

Thirty  material  properties  are  required.  Ten  of  these  properties  are  global  and  the 
remaining  twenty  are  associated  with  each  of  the  four  mechanisms.  The  global  properties  are 
listed  in  Table  4.2  and  the  mechanism-specific  properties  are  listed  in  Table  4.3. 

Associated  Parameters  -  Global  Parameters 

Two  pairs  of  global  parameters  are  associated.  One  pair  is  used  to  adjust  the  friction 
angle  for  the  effects  of  mean  stress.  Then,  a  yield  criterion  is  determined  based  on  adjusted 
friction  angle.  The  parameters  are  phi  ratio  (PHIRATIO)  and  an  Over  Consolidation  factor 
(Decay).  The  expression  used  in  the  code  for  OC  factor  is  the  ratio  of  the  reference  stress,  Pe, 
to  the  mean  stress,  am. 

The  expression  for  the  Yield  limit  is  based  on  the  formulation  11*12/13,  where  the  “I” 
terms  are  the  stress  invariants.  The  stress  tensor  used  in  the  calculation  is  given  in  vector  form 
as  STRESS  (6).  Stress  1  has  the  magnitude  of  (1  +  sin  (<(>))/(l-sin  (<)>)).  The  shearing  stresses 
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are  zero  (STRESS  (4),  STRESS  (5),  and  STRESS  (6)).  The  normal  stresses  are  given  as 
principal  stresses.  STRESS  (2  )and  STRESS  (3)  are  unity. 

The  other  pair  of  associated  parameters  is  used  in  the  shear-volume  coupling  term. 
The  volume  change  is  proportional  to  the  plastic  strain.  The  volume  change  is  the  difference 
of  two  terms.  The  first  term  is  the  inner  product  of  the  shear  stress  and  the  total  plastic  strain 
and  that  quantity  normalized  by  the  mean  stress.  The  second  term  is  square  root  of  inner 
product  of  the  plastic  strains.  This  quantity  is  a  scalar  that  is  the  magnitude  of  the  plastic 
strain.  This  term  is  multiplied  by  a  dilatancy  factor,  Me.  In  this  model  a  scaling  factor, 
gamma,  was  introduced  to  reduce  the  effect  of  the  shear  volume  coupling.  Gamma  is  unity  in 
traditional  CSSM. 

Associated  Parameters  -  Mechanism  Parameters 

Each  mechanism  acts  without  consideration  of  the  other  mechanisms.  That  is, 
subroutine  Ammos  is  called  once  per  mechanism  and  performs  its  calculations  without 
consideration  of  previous  calls.  However,  selection  of  mechanism  parameters  does  require 
some  consideration  of  all  of  mechanisms  acting  as  a  unit  of  four.  The  stiffness,  both  shear  and 
volumetric,  must  be  distributed  among  the  mechanisms  such  that  the  sum  of  each  mechanism 
stiffness  equals  the  global  stiffness  parameters. 

Flow  Scheme 

The  main  UMAT  initializes  variables  and  calls  Subroutine  Sand_driver  seven  times. 
The  first  call  to  Sand_driver  returns  a  solution  for  the  stress  and  updated  solution-dependent 
state  variables  (SDV’s).  The  remaining  calls  return  portions  of  the  Jacobian.  The  main 
UMAT  updates  the  Jacobian,  stress,  and  SDV’s. 

Solution  Dependent  Variables  (SDV’s) 

The  constitutive  model  has  internal  variables  whose  purpose  is  similar  to  the  internal 
variables  of  endochronic  theory.  The  internal  variables  hold  the  stress  state  of  each 
mechanism  in  terms  of  an  internal  force.  Two  types  of  internal  forces  are  used.  Qs  are  the 
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shear  forces,  and  Qh  are  the  hydrostatic  forces.  Each  mechanism  has  six  shear  forces  and  one 
hydrostatic  force  associated  with  its  strain  history.  A  total  of  seven  internal  forces  are  needed 
for  each  mechanism.  The  four  mechanisms  require  28  internal  forces  to  be  carried  through 
each  step.  The  void  ratio  must  be  carried  as  well.  Twenty-nine  SDV’s  are  used  in  all. 

UMAT  Main  Subroutine 

The  UMAT  is  called  from  the  ABAQUS  program,  herein  described  as  the  main 
calling-program.  The  UMAT  main  program  initializes  stress  to  values  sent  from  the 
ABAQUS  main  calling-program  and  assigns  properties  to  values  set  in  the  UMAT  control 
card.  The  internal  variables  (SDV’s),  total  strain  and  strain  increments  enter  the  UMAT  with 
the  values  passed  by  the  main  calling-program.  A  flow  chart  for  the  UMAT  main  program  is 
shown  in  Figure  4.3. 

Then,  the  UMAT  main  program  calls  subroutine  Sand_driver  passing  all  of  the 
stresses,  strains,  strain  increments,  and  internal  variables.  Sand_driver  returns  appropriate 
stress  and  internal  variables  that  may  have  been  clipped  if  the  material  yielded.  The  main 
subroutine  pushes  the  new  stresses  and  internal  variables  into  the  appropriate  arrays  and  then 
prepares  dummy  strains  to  send  to  Sand_driver  for  the  purpose  of  determining  the  Jacobian. 
Sand_driver  is  called  again  in  a  loop  to  create  the  data  for  the  six  Jacobian  terms.  Once  the 
loop  is  complete,  the  main  program  returns  the  Jacobian,  the  updated  stresses  and  updated 
internal  variables  to  the  ABAQUS  main  calling-program. 

Subroutine  Sand_driver 

Sand_driver  is  called  from  the  UMAT  main  program.  Sand_driver  is  provided  with 
strains,  and  internal  variables,  and  stresses.  Sand_driver  calculates  plastic  strain,  volumetric 
strain,  computes  a  normalizing  stress  from  the  NCL  variable  and  the  void  ratio,  and 
determines  a  hydrostatic  parameter  associated  with  the  internal  variables  for  hydrostatic  stress. 
Also,  the  internal  variables  for  each  mechanism  are  updated.  Flow  charts  for  the  Sand_driver 
subroutine  are  shown  in  Figures  4.4  and  4.5.  The  yield  limit  is  calculated  for  each 
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mechanism.  The  yield  limit  is  a  function  of  the  frictional  strength  of  each  mechanism. 
Sand_driver  calls  two  other  subroutines,  Ammos  and  Hydros  in  order  to  perform  these 
calculations.  Sand_driver  returns  to  the  UMAT  main  with  updated  stresses  and  internal 
variables. 

Subroutine  Sand_driver  is  called  again  by  a  loop  in  the  UMAT  main  to  determine  data 
for  the  tangent  stiffness  matrix  or  Jacobian.  The  UMAT  main  provides  updated  stress  and 
internal  variables  and  dummy  strain  to  enable  a  partial  derivative  to  be  estimated  for  each  term 
of  the  tangent  stiffness  matrix  or  Jacobian. 

Subroutine  Ammos 

Ammos  is  used  to  set  a  yield  limit  and  check  the  shearing  stress  produced  by  the 
incremental  strains  provided  in  the  call  from  the  Sand_driver  subroutine.  Ammos  is  called 
from  Sand_driver  one  time  for  each  of  the  four  mechanisms.  Ammos  is  sent  incremental 
strains,  internal  variables  and  yield  limit  for  the  mechanism.  Ammos  checks  for  a  mean 
tensile  stress  and  sets  the  value  of  mean  stress  to  a  small  compressive  value  if  tension  was 
detected.  The  shear  strains  are  determined  and  the  shear  stress  increment  is  determined 
assuming  that  the  strain  was  elastic.  The  location  of  the  yield  surface  is  determined  for  the 
mechanism  based  on  the  values  of  the  internal  variables  and  compared  to  the  yield  limit.  A 
clipping  subroutine  RadialRetum  is  called  if  the  shear  stress  point  is  located  beyond  the  yield 
limit.  If  clipping  was  required  due  to  yielding,  Ammos  updates  the  plastic  strain  for  each 
mechanism  and  the  total  plastic  strain.  Ammos  records  the  plastic  strain  as  zero  if  no  clipping 
was  necessary.  Ammos  updates  and  returns  the  values  of  shear  stress  and  shear  internal 
variables  to  Sand_driver.  A  flow  chaft  for  the  Ammos  subroutine  is  shown  in  Figure  4.6. 
Subroutine  Hydros 

Hydros  is  called  by  Sand_driver  and  used  to  update  the  hydrostatic  internal  variables 
and  clip  the  hydrostatic  internal  variables  if  either  the  compression  limit  or  the  tensile  limit 
were  exceeded.  Since  the  hydrostatic  stress  can  be  described  as  a  scalar  quantity,  this 
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subroutine  is  much  simpler  than  AMMOS.  Hydros  returns  the  updated  hydrostatic  internal 
variables  (hydrostatic  stresses)  to  Sand_driver. 

Subroutine  RadialReturn 

RadialRetum  is  a  clipping  subroutine  called  by  Amnios.  The  subroutine  performs 
radial  return  of  stress  point  to  yield  function  (Matsuoka  and  Nakai,  1977),  given  by: 

Fy(Q)  =  11*12/13  (4.5) 

where,  II,  12,  and  13  are  the  stress  invariants.  A  transformation  is  first  performed  to  principal 
stress  space,  then  the  return  is  performed  such  that  II  and  (Pv2-Pv3)/(Pvl-Pv3)  are  held 
constant.  Pvl,  Pv2,  and  Pv3  are  the  principal  stress  values.  With  these  constraints,  Fy  = 
Ylimit  becomes  a  cubic  equation.  The  stress  tensor  is  computed  from  the  eigenvectors  and 
adjusted  eigenvalues.  Therefore,  the  adjusted  stress  tensor  has  the  same  principal  axes,  mean 
stress,  and  Lode  parameter  as  the  original  stress  tensor. 

Summary  of  Calling  Schedule 

For  each  time  the  UMAT  is  called.  Subroutine  Sand_driver  is  called  seven  times, 
once  for  the  stresses  and  internal  variables  and  six  times  for  the  Jacobian  or  tangent  stiffness 
matrix  required  by  ABAQUS.  Subroutine  Ammos  is  called  eight  times  per  call  to  Subroutine 
Sand_driver.  Subroutine  Hydros  is  called  four  times  per  call  to  Sand_driver.  Subroutine 
RadialRetum  could  be  called  a  maximum  of  one  time  per  call  to  Ammos.  Subroutine 
RadialRetum  is  called  only  when  plastic  strain  has  occurred.  Table  4.4  shows  the  range  of 
potential  numbers  of  calls  to  each  subroutine  per  iteration  of  each  load  increment  of  each  step 
in  an  ABAQUS  analysis  (Meade,  1998). 

Model  Operation 

In  order  to  demonstrate  the  operational  characteristics  of  the  WES  Multimechanical 
constitutive  model  (WES  MM)  the  following  discussion  of  a  cyclic  stress  strain  curve  is 
presented.  An  idealized  representation  of  the  WES  MM  model  is  shown  in  Figure  4.6.  The 
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major  points  emphasized  are  that  the  model  has  four  mechanisms  with  elastic-plastic  behavior. 
The  strain  experienced  under  load  application  is  common  to  all  elements,  and  the  stiffness  and 
yield  of  each  element  are  different.  Figure  4.7  shows  the  stress  strain  curve  for  a  cyclic  test 
with  the  crucial  stages  of  the  test  numbered  as  Points  1-11. 

The  stress  path  from  which  the  above  stress-strain  curve  is  derived  is  presented  in 
Figures  4.8  through  4. 1 8.  A  separate  figure  for  each  of  the  critical  points  shown  in  Figure  4.7 
is  used  to  describe  the  various  stages  of  yielding  and  stress  reversals.  In  many  plasticity 
models,  a  hardening  law  is  employed  to  describe  the  change  in  yield  strength  that  accompanies 
the  occurrence  of  plastic  strain.  From  the  discussion  presented  in  Chapter  3  it  is  evident  that 
these  hardening  rules  can  become  very  complex  and  difficult  to  implement.  The  WES  MM 
model  employs  4  predefined  yield  surfaces  to  capture  the  hardening  that  occurs  in  a  material 
loaded  beyond  an  initial  yield  stress. 

At  point  1  the  first  mechanism  has  yielded  and  begun  experiencing  plastic 
deformation.  Figure  4.8  shows  the  stress  path  of  each  of  the  four  mechanisms  at  Point  1  in  a 
principal  stress  difference  versus  mean  normal  stress  space  ( q  versus  p).  The  stress  path  of 
each  of  the  mechanisms  is  shown  in  a  separate  plot  labeled  Ml  through  M4.  The  yield 
surfaces  for  the  mechanisms  are  shown  as  the  thin  lines  emanating  from  the  origin  of  each  q-p 
axis.  The  actual  stress  path  of  the  elements  is  shown  as  the  dark  lines  moving  off  the 
horizontal  axis  at  some  distance  p  from  the  origin.  The  stress  difference  or  shear  stress  in  a 
yielded  mechanism  increases  after  yield  only  as  a  function  of  the  increase  of  normal  stress,/?. 

At  Point  2  Mechanism  2  has  yielded  and  begun  accumulating  plastic  strain  along  with 
Mechanism  1 .  The  stress  path  at  Point  2  is  shown  in  Figure  4.9.  Mechanisms  3  and  4 
continue  to  respond  elastically  until  Mechanism  3  yields  at  Point  3  (Figure  4.10). 

At  Point  4  (Figure  4.1 1)  the  stress  is  reversed  and  unloading  begins.  During  the  initial 
stage  of  unloading,  all  mechanisms  are  undergoing  elastic  strain.  As  the  mechanisms  reach 
yield  in  extension  (Figures  4.12  and  4.13),  the  stress-strain  curve  breaks  over  to  change  slope 


49 


just  as  it  does  in  loading.  The  mean  stress  is  higher  at  each  of  the  breaks  points  in  this  unload- 
reload  cycle  and  therefore  the  shear  stress  at  yield  is  higher. 

Figure  4.15  depicts  the  stress  path  at  the  time  of  yielding  in  Mechanism  1  under 
reloading.  In  Figure  4.7  this  occurs  at  Point  8  with  a  higher  yield  stress  than  that  seen  in  the 
initial  loading  curve.  The  same  type  of  behavior  is  seen  in  Mechanism  2  at  Point  9  as  shown 
in  Figure  4.16. 

At  Point  9  in  Figure  4.7,  the  third  mechanism  yields  in  reload  as  shown  in  Figure  4. 17. 
Again  this  occurs  at  a  higher  mean  and  shear  stress  than  the  initial  loading.  The  resulting 
hysteresis  loops  formed  from  Points  4-10  produce  permanent  deformation  under  cyclic 
loading  conditions.  Figure  4.18  shows  the  continued  loading  resulting  in  plastic  strains  for  all 
but  Mechanism  4,  which  remains  elastic.  This  ratcheting  effect  produces  a  strain  that 
increases  with  load  repetitions. 
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Table  4.2.  Global  Properties 


Name 

Label  in  code 

Comments 

Phi 

PHILIMIT 

friction  angle 

Cohesion 

C 

cohesion 

Bulk  Modulus 

K 

Shear  Modulus 

G 

phi  ratio 

PHIRATIO 

Hydrostatic  Intercept 

Fh 

Intercept  of  Normal 

Consolidation  Line  (NCL) 

Reciprocal  of  Cc 

BETA 

Reciprocal  of  the  slope  of  NCL 

Shear-volume  factor 

Me 

shear-volume  coupling  term 

OC  factor 

Decay 

strength  reduction  term 

dilatancy  scaling  factor 

GAMMA 

Table  4.3.  Mechanism  Properties 


Name 

Label  in  code 

Comments 

Strength  factor 

PHIFRAC 

scales  friction  angle 

Mean  Stress  factor 

PFACT 

scales  mean  stress 

Shear  Stiffness  factor 

SHEARRATIO 

distributes  shear  stiffness 

Compression  limit 

HLIMIT 

absolute  compression  limit 

Volumetric  Stiffness  factor 

BULKRATIO 

distributes  volumetric  stiffness 

Table  4.4.  Frequency  of  Calls 


Subroutine 

Relative  #  of  calls 

Total  #  of  calls  per  Step 

Sand_driver 

7  per  call  to  UMAT 

7 

Ammos 

8  per  call  to  Sand_driver 

56 

RadialRetum 

0  or  1  per  call  to  Ammos 

0  to  56 

Hydros 

4  per  call  to  Sand_driver 

28 

Subroutine 
SDVINI 
initializes 
state  variable 
array  on  first 
pass  through 
UMAT 


Strain  Increment  Tensor 
Material  Properties 
State  Variable  Array 


/ 

/ 

1 

r 

Call  Sand_Driver  to  Calculate 

Updated  Stress  Tensor 
Plastic  Strain  Tensor 

r 

Loop  6  times 


Exit 


Call  Sand  Driver  to  Calculate  Jacobian  for  ABAQUS 

. 

Return  to  ABAQUS  with  Updated  Stress 
Tensor,  State  Variable  Array,  and  Jacobian 


Figure  4.3.  Flow  chart  for  WES  MM  ABAQUS  UMAT 


52 


53 


Continued 


Figure  4.5.  Flow  chart  for  Subroutine  Sand_driver  (Part  2) 
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AMMOS  Subroutine 

r  - \ 

Subroutine  AMMOS 


Return  to  Sand_driver  with  Updated 
Plastic  Strains  and  Shear  Stress  Array 


Figure  4.6.  Flow  chart  for  Subroutine  Ammos 
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Common  Strain 


Figure  4.8.  Stress  versus  strain  for  a  cyclic  test 
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CHAPTER  5:  MODEL  CALIBRATION 


GENERAL 

In  order  to  properly  apply  any  constitutive  model  to  predict  the  response  of  materials 
under  load  the  models  must  be  calibrated  with  test  data.  In  essence  the  parameters  used  to 
define  strength,  failure,  and  deformation  properties  must  be  defined  for  any  material  to  be 
modeled.  This  chapter  describes  the  model  requirements,  laboratory  tests,  and  analysis  to 
achieve  a  proper  calibration  for  both  the  ABAQUS  Drucker-Prager  model  and  the  WES 
Multimechanical  model. 

ABAQUS  DRUCKER-PRAGER  MODEL 

The  model  uses  three  groups  of  parameters:  stiffness  parameters,  failure  surface 
parameters,  and  cap  parameters.  The  general  procedures  used  to  determine  these  parameters 
from  laboratory  test  data  are  presented  in  this  section. 

Failure  Surface 

As  presented  in  Chapter  4,  the  ABAQUS  Drucker-Prager  failure  surface  is  written  in  a 
q  (principal  stress  difference)  versus p  (mean  normal  stress)  space  as: 

Fs  =  q  - ptan  J3 -d  =  0  (4.1) 

where  /?  and  d  represent  the  angle  of  friction  of  the  material  and  its  cohesion,  respectively. 

Cap  Yield  Surface 

The  cap  yield  surface  has  an  elliptical  shape  with  constant  eccentricity  in  q-p  plane 
and  also  includes  dependence  on  the  third  stress  invariant  in  the  deviatoric  plane.  The  cap 
surface  hardens  or  softens  as  a  function  of  the  volumetric  inelastic  strain.  The  ABAQUS 
Drucker-Prager  failure  surface  is  written  in  a  q  (principal  stress  difference)  versus  p  (mean 
normal  stress)  space  as  Equation  4.2, 

Fc  =J[P-Paf  + 


Rq 


(1  +  a -a /cos  (3) 


-R(d  +  pa  tan/?)  =  0 


(4.2) 
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where  R  is  a  material  parameter  that  controls  the  shape  of  the  cap,  a,  a  cap  transition  factor), 
and  pa  is  an  evolution  parameter  that  represents  the  volumetric  inelastic  strain  driven 
hardening/softening.  The  pa  parameter  is  a  function  of  the  plastic  volumetric  strain  and 
volumetric  yield  stress.  The  materials  typically  used  in  granular  base  courses  in  pavements 
have  a  very  high  level  of  compaction  and  strength.  One  would  only  expect  to  intersect  the  cap 
in  such  materials  under  loads  much  higher  than  those  experienced  in  pavements,  such  as  blast 
or  shock  conditions.  In  essence  this  reduces  the  cap  model’s  operation  back  to  a  simpler  two- 
parameter  friction  model  based  on  /?  and  d.  Figure  5.1  shows  the  simplified  mode  parameters 
with  the  stress  regime  of  interest  in  the  shaded  area. 


Transition 
surface,  Ft 


Figure  5.1.  ABAQUS  Drucker-Prager  model  with  stress  regime  of  interest  shown  in  gray 
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Calibration 


At  least  three  experiments  are  required  to  calibrate  the  simplest  version  of  the  Drucker 
-Prager  model:  a  hydrostatic  compression  test  and  two  triaxial  compression  tests  (more  than 
two  tests  are  useful  for  a  more  accurate  calibration). 

The  hydrostatic  compression  test  is  performed  by  pressurizing  the  sample  equally  in 
all  directions.  The  applied  pressure  and  the  volume  change  are  recorded.  Triaxial 
compression  experiments  are  performed  using  a  standard  triaxial  machine  where  a  fixed  con¬ 
fining  pressure  is  maintained  while  the  differential  stress  is  applied.  Several  tests  covering  the 
range  of  confining  pressures  of  interest  are  usually  performed.  Again,  the  stress  and  strain  in 
the  direction  of  loading  are  recorded,  together  with  the  lateral  strain  so  that  the  correct  volume 
changes  can  be  calibrated.  Unloading  measurements  in  these  tests  are  useful  in  determining 
elastic  properties,  particularly  in  cases  where  the  initial  elastic  region  is  not  well  defined. 

The  stress-strain  curve  from  the  hydrostatic  compression  test  gives  the  evolution  of 
the  hydrostatic  compression  yield  stress.  The  friction  angle,  /?,  and  cohesion,  d,  which  define 
the  shear  failure  dependence  on  hydrostatic  pressure,  are  calculated  by  plotting  the  failure 
stresses  of  any  two  uniaxial  and/or  triaxial  compression  experiments  in  q  (principal  stress 
difference)  versus  p  (mean  normal  stress)  space:  the  slope  of  the  straight  line  passing  through 
the  two  points  gives  the  angle  /?  and  the  intersection  with  the  9-axis  gives  d. 

WES  MULTIMECHANICAL  CONSTITUTIVE  MODEL 
The  model  uses  three  groups  of  parameters:  stiffness  parameters,  strength  parameters, 
and  a  shear-volume  coupling  parameter.  The  general  procedures  used  to  determine  these 
parameters  from  laboratory  test  data  are  presented  in  this  section. 

Calibrating  the  Model  -  General  Approach 
The  procedure  for  calibrating  the  model  requires  a  set  of  several  triaxial  tests,  either 
drained  or  undrained  with  pore  pressure  measurements.  First,  the  relation  between  mean 
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effective  stress  and  void  ratio  at  the  ultimate  state  is  plotted  similar  to  an  e  -  log  p  curve  as 
shown  in  Figure  5.2.  The  slope  Cc  and  intercept,  Fh,  are  used  to  determine  the  relation 
between  void  ratio  and  the  reference  pressure,  Pe.  Next,  the  hydrostatic  stress-strain  curve  is 
plotted  in  a  normalized  form  in  which  the  hydrostatic  stress  is  divided  by  the  reference  stress. 
In  this  form  the  hardening  effect  of  void  ratio  decrease  is  removed,  leaving  the  fundamental 
curve.  The  normalized  curve  is  then  divided  into  regions  to  be  represented  by  each 
mechanism.  The  yield  stress  associated  with  each  mechanism  is  thus  determined.  The  stiffness 
of  each  mechanism  is  determined  by  the  change  in  modulus  that  occurs  as  each  yield  limit  is 
crossed. 

A  similar  procedure  is  carried  out  for  the  shear  response.  The  shear  yield  limit  is 
determined  for  each  mechanism.  Friction  angles  are  selected  based  on  the  ultimate  friction 
angle  at  a  stress  level  close  to  that  of  the  expected  service  loads.  From  these  data,  the 
distribution  factor  for  hydrostatic  stress  can  be  determined  for  each  mechanism.  The 
calibration  from  the  shear  moduli  is  the  same  as  that  for  the  bulk  moduli  of  the  hydrostatic 
mechanism. 


Figure  5.2.  Void  ratio  versus  log  normal  stress  plot  used  to  determine  NCL  for  WES  MM 
model 
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LABORATORY  TESTS 


Tests  were  conducted  on  a  well-graded  limestone  base  course  material  to  determine  its 
response  to  loads  and  to  define  its  yield  surface  for  use  with  plasticity  formulations  such  as  the 
Mohr-Coulomb  and  Drucker-Prager  models.  Five  Unconfined  Compression  (UCC)  tests  were 
conducted.  Conventional  Triaxial  Compression  (CTC)  tests  were  conducted  at  four  confining 
pressures  up  to  80  psi  (551.6  kPa),  with  axial  strains  up  to  5  percent.  Uniaxial  Strain  (UXE) 
tests  and  Hydrostatic  Compression  (HC)  tests  were  conducted  up  to  confining  pressure  levels 
of  100  psi  (689.5  kPa).  Replicates  of  each  test  were  performed  to  insure  that  variations  in 
response  could  be  identified  and  corrected. 

The  mechanical  response  of  granular  materials  must  be  clearly  understood  to 
accurately  predict  the  performance  of  flexible  pavements.  Due  in  large  part  to  testing 
difficulties,  the  measurement  of  load-induced  response  of  granular  materials  has  received  little 
attention  in  the  geotechnical.  Considerable  effort  expended  during  this  research  was  aimed  at 
developing  equipment,  procedures  and  skills  necessary  for  preparing  and  testing  unbound, 
highly  angular,  granular  materials. 

Material 

The  Type  610  (MDOT,  1990)  well-graded  crushed  limestone  material,  as  shown  in 
Figure  5.3,  selected  for  this  study  was  used  as  a  base  course  in  an  airfield  pavement  test 
section  (Webster  1993).  The  grain-size  analysis  is  shown  in  Figure  5.4.  This  grain  size 
distribution  is  typical  for  aggregate  base  course  materials  used  in  many  airfield  pavements. 
Webster  reported  the  material  to  be  an  SW-SC  according  to  the  Unified  Soil  Classification 
System.  However,  further  investigation  proved  that  the  material  was  actually  a  GW  material. 

It  had  a  liquid  limit  of  17,  a  plastic  limit  of  1 1,  and  a  plasticity  index  of  6.  Using  modified 
proctor  procedures,  in  accordance  with  ASTM  D  1557,  optimum  moisture  content  for 
compaction  was  determined  to  be  4.5  percent.  Dry  unit  weight  at  optimum  moisture  content 
was  determined  to  be  144  lb/ft3  (2306.7  kg/m3). 
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Prior  to  performing  each  mechanical  property  test,  the  height,  diameter  and  weight  of 
each  remolded  specimen  were  determined.  These  measurements,  along  with  the  aggregate 
specific  gravity  and  water  content,  were  used  to  calculate  dry  density  and  void  ratio  for  each 
specimen.  The  variation  of  height,  weight,  and  diameter  of  the  specimens  were  carefully 
controlled  to  arrive  at  the  dimensions  shown  in  Table  5.1.  All  specimens  were  constructed 
using  the  GTM  procedure  to  the  dimensions  and  the  weight  shown.  Specimens  were  fabricated 
to  reproduce  the  field  density  of  137.2  lb/  ft3  (2199.5  kg/m3)  and  moisture  content  of  4.0  %. 

A  digital  electronic  caliper  with  accuracy  of  +/-  0.001  inches  (0.03  mm)  was  used  to  verify 
specimen  dimensions  prior  to  testing.  A  digital  electronic  scale  with  a  maximum  range  of  22 
lbs.  (10  kg)  and  an  accuracy  of  +/-  0.0002  lbs.  (0.1  g)  was  used  to  verify  specimen  weight 
prior  to  testing.  Specimens  not  meeting  weight  and  dimension  requirements  were  rejected  for 
testing. 


Table  5.1.  Granular  Limestone  Specimen  Properties 


Diameter 
In.  (mm) 

Height 
in.  (mm) 

Area 

in2 

(mm2) 

Dry 
Weight 
lb.  (kg) 

Volume 

ft3 

(m3) 

Density 

lb/ft3 

(kg/m3) 

Void 

Ratio 

Moisture 

Content 

% 

4.00 

(101.6) 

8.36 

(213.4) 

8.38 

(3.81) 

0.061 

(0.00173) 

137.2 

(2199.5) 

0.21 

4.0 

Specimen  Preparation 

Each  granular  specimen  tested  in  this  study  was  exactly  4  inches  (101.6  mm)  in 
diameter  by  8.36  inches  (213.4  mm)  in  height  as  seen  in  table  5.1.  The  compacted  specimens 
were  tested  in  a  conventional  triaxial  compression  chamber  meeting  ASTM  D2850,  which 
consisted  of  a  reinforced  Plexiglas  pressure  vessel,  a  stainless  steel  base,  and  a  stainless  steel 
top.  All  specimens  were  tested  in  the  triaxial  chamber  immediately  following  completion  of 
the  compaction  process  to  insure  that  no  damage  or  moisture/strength  loss  occurred  during 
extended  storage  periods. 


73 


Percent  Finer  by  Weight 


Figure  5.3.  Well-graded  crushed  limestone  used  in  laboratory  tests 


Sieve  Opening,  inches  Sieve  Numbers 


Grain  Size  in  Millimeters 


Figure  5.4.  Grain-size  analysis  of  well-graded  crushed  limestone  used  in  laboratory  tests 
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Specimens  were  prepared  using  the  Corps  of  Engineers  Gyratory  Compaction  Testing 
Machine  (GTM)  (ASTM  D  3387).  Compaction  of  materials  using  the  gyratory  method 
applies  normal  forces  to  both  the  top  and  bottom  faces  of  the  material  confined  in  cylindrically 
shaped  molds.  Normal  forces  at  designated  pressures  are  supplemented  with  a  kneading 
action  or  gyratory  motion  to  compact  the  material  into  a  denser  configuration  with  aggregate 
particle  orientation  more  consistent  with  in-place  pavements. 

The  gyratory  compaction  method  involves  placing  loose  material  into  a  4-inch 
diameter  by  10-inch  length  mold  and  loading  into  the  GTM  at  a  prescribed  normal  stress  level 
which  represents  anticipated  traffic  contact  pressure.  The  material  and  mold  are  then  rotated 
through  a  1 -degree  gyration  angle  for  a  specified  number  of  revolutions  of  the  roller  assembly. 
This  compaction  process  produces  stress-strain  properties  that  are  representative  of  those  in  a 
field  compacted  material  (Ahlrich,  1997).  A  schematic  of  the  gyratory  compaction  device  is 
shown  in  Figure  5.5. 

The  gyratory  testing  machine  shown  in  Figure  5.6  was  used  to  compact  all  laboratory 
specimens  in  this  research.  The  gyratory  compactive  effort  used  in  this  laboratory  study  was  a 
200-psi  (1378.7  kPa)  normal  stress  level,  1-degree  gyration  angle  and  30-  50  revolutions  of 
the  roller  assembly.  This  compaction  effort  produced  specimens  that  were  nominally  8  inches 
long  and  4  inches  in  diameter.  The  specimens  were  sealed  with  1-in.  (25.4-mm)  thick 
aluminum  endcaps  and  double  0.025-in.  (0.635-mm)  latex  membranes  before  being  place  in 
the  triaxial  testing  device.  A  target  density  of  137.2  lb/ft3  (2199.5  kg/mm3 )  was  used  to  select 
the  compaction  effort  described  above.  The  gyratory  compaction  process  produced  highly 
repeatable  samples  and  contributed  greatly  to  the  success  of  the  laboratory  testing  phase  of 
this  research.  Experiences  with  other  compaction  methods  such  as  vibratory  and  hammer 
compaction  procedures  for  granular  materials  proved  unsuitable  for  this  investigation. 
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Figure  5.5.  Schematic  of  gyratory  testing  machine 


Figure  5.6.  Gyratory  testing  machine  used  for  specimen  preparation 
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Description  of  Test  Device 


A  conventional  cylindrical  soils  triaxial  testing  device  conforming  to  ASTM  D2850 
was  used  to  perform  the  mechanical  property  tests.  The  test  device  had  overall  nominal 
dimensions  18  inches  (457.2  mm)  in  height  by  12  inches  (304.8  mm)  in  diameter  with  a 
capacity  to  test  specimens  up  to  5  inches  (12.7  mm)  in  diameter.  The  pressure  vessel  was 
reinforced  Plexiglas  with  hardened  stainless  steel  encaps  and  connecting  rods. 

The  confining  pressure  was  supplied  by  air  pressure.  A  servo-controlled  Instron 
testing  machine,  capable  of  applying  tensile  or  compressive  loads  up  to  60,000  lb.  (266  kN) 
supplied  the  axial  load.  The  Instron  testing  machine  and  triaxial  chamber  is  shown  in  Figure 
5.7.  The  loader  could  be  controlled  either  manually  or  by  computer  in  order  to  produce  a 
desired  rate  of  loading  or  displacement.  The  input  to  the  servo-control  unit  was  produced  by  a 
function  generator,  which  could  be  programmed  to  produce  large  variety  of  load  or 
displacement  histories.  A  load  cell  measured  the  axial  force  applied  to  each  granular 
specimen.  The  confining  pressure  applied  to  the  specimens  is  measured  with  a  pressure 
transducer,  located  at  the  air  supply  regulator. 

Measurement  of  the  changes  in  the  specimen  dimensions  were  critical  considerations 
in  the  testing.  Measurements  of  deformation  under  load  of  a  remolded  cohesionless  material 
is  a  very  difficult  task.  The  measurement  devices  must  provide  for  accurate  changes  in  length 
and  diameter  without  affecting  the  response  of  the  material. 

Changes  in  specimen  length  were  measured  with  two  diametrically  opposed  linear 
variable  differential  transducers  (LVDT),  mounted  on  the  end  platens  on  the  inside  of  the 
chamber.  The  change  in  diameter  of  the  specimen  under  load  was  measured  with  a  device  that 
consisted  of  four  strain-gaged  spring  arms  attached  to  a  mounting  ring  and  calibrated  to 
provide  a  diameter  change  output  in  a  full  bridge  configuration.  A  photo  of  a  specimen  with 
its  deformation  devices  attached  is  shown  in  Figure  5.8. 
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Figure  5.7.  Instron  servo-controlled  testing  machine 
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Figure  5.8.  Granular  limestone  specimen  with  instrumentation  attached 
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RESULTS  OF  LABORATORY  TESTS 


Unconfined  Compression  Tests 

Five  unconfined  compression  tests  were  conducted  on  remolded  4-in  by  8-in 
specimens.  The  tests  were  conducted  in  the  same  chamber  as  the  CTC  tests.  Each  unconfined 
compression  (UC)  test  was  conducted  by  applying  an  axial  load  with  a  constant  rate  of  1%  per 
minute.  The  load  was  applied  until  the  granular  material  exhibited  either  a  maximum  axial 
stress  (peak)  or  an  axial  strain  of  5  percent.  With  a  specimen  length  of  8.36  inches  (212.3 
mm)  a  strain  of  5  percent  equals  a  total  change  in  length  of  0.418  inches  (10.  9  mm). 
Controlling  the  test  by  a  specified  rate  of  strain  allowed  for  the  capture  of  softening  response 
during  post-yield  stress  application. 

A  summary  of  these  tests  is  shown  in  Table  5.2.  Plots  of  Axial  stress  versus  axial 
strain  are  shown  in  Figure  5.9.  Figure  5.10  shows  a  plot  of  mean  normal  stress  versus 
volumetric  strain  for  the  tests.  The  individual  plots  of  test  data  from  the  unconfined 
compression  tests  are  shown  in  Appendix  E. 


Table  5.2.  Summary  of  Results  at  Maximum  Axial  Stress  from  UCC  Tests 


Axial 

Strain 

Radial 

Strain 

Volumetric 

Strain 

Axial  Stress 

Mean  Normal 
Stress 

% 

% 

% 

psi 

kPa 

psi 

kPa 

Ucc_l 

1.9 

-6.0 

-10.1 

m 

52.7 

2.5 

17.6 

Ucc_2 

1.3 

-4.2 

-7.1 

8.4 

57.6 

2.8 

19.2 

Ucc_3 

1.3 

-6.0 

11.7 

80.7 

3.9 

26.9 

Ucc_4 

1.9 

-5.2 

-8.5 

10.1 

69.6 

23.2 

Ucc_5 

2.3 

-6.1 

-9.8 

8.8 

61.0 

2.9 

20.3 

Mean 

1.7 

-5.0 

-8.3 

9.3 

64.3 

3.1 

21.4 

0.4 

1.1 

1.8 

1.6 

11.0 

0.5 

*  Standard  Deviation 
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Crushed  Limestone  Typa  610 
Unconfined  Compression 


Axial  Strain,  % 


Figure  5.9.  Axial  stress  versus  strain  for  unconfined  compression  tests  of  granular  limestone 
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Figure  5.10.  Mean  normal  stress  versus  volumetric  strain  for  unconfined  compression  tests  of 
granular  limestone 
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Conventional  Triaxial  Compression  Tests 


Conventional  triaxial  tests  (CTC)  were  conducted  according  to  ASTM  D2850  except 
for  the  displacement  measuring  system  and  subtle  differences  required  to  prevent  damage  of 
the  specimen  dining  assembly  of  the  device.  Tests  at  four  confining  pressure  levels  with  a  at 
least  three  repetitions  at  each  level  of  confining  pressure  were  conducted.  Each  conventional 
triaxial  compression  (CTC)  test  was  conducted  in  two  phases.  An  isotropic  compression  (IC) 
phase  was  conducted  by  applying  a  confining  pressure  to  all  sides  of  the  cylindrical  specimen, 
while  measuring  its  change  in  height  and  diameter.  These  data  are  often  plotted  as  mean 
normal  stress  versus  volumetric  strain,  the  slope  of  which  is  the  bulk  modulus,  K.  After  the 
desired  confining  pressure  had  been  attained  during  the  IC  phase,  the  triaxial  compression 
phase  was  conducted.  This  was  accomplished  by  applying  an  axial  load  with  a  strain  rate  of 
1%  per  minute,  while  the  confining  pressure  was  held  constant.  After  the  maximum  strain  of 
5  percent  was  reached  the  test  machine  was  reversed  to  allow  measurement  of  unloading 
response.  These  tests  were  essentially  undrained  tests  that  did  not  generate  any  excess  pore 
water  pressure.  The  specimens  had  saturation  levels  in  the  range  of  50%  with  void  ratios  on 
the  order  of  0.21  at  the  beginning  of  each  test. 

A  summary  of  these  tests  is  shown  in  Table  5.3.  The  data  in  Table  is  5.3  is  organized 
according  to  confining  pressure  with  a  statistical  summary  of  each  level  of  response  provided 
in  the  table.  Plots  of  principal  stress  difference  versus  principal  strain  difference  are  shown  in 
Figure  5.11.  Plots  of  principal  stress  difference  versus  mean  normal  stress  are  shown  in 
Figure  5.12.  Figure  5.13  shows  a  composite  plot  of  mean  normal  stress  versus  volumetric 
strain.  The  individual  plots  of  test  data  from  the  conventional  triaxial  compression  tests  are 
shown  in  Appendix  C.  The  tests  are  designated  as  CTCxx_y.  The  coding  designation  xx  is 
the  confining  pressure  in  pounds  per  square  inch,  and  y  is  the  replicate  number  at  the  confining 
pressure  xx.. 
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Table  5.3.  Summary  of  Results  at  Maximum  Axial  Stress  from  CTC  Tests 


Axial 

Strain 

Radial 

Strain 

Strain 

Difference 

Confining 

Pressure 

Volumetric 

Strain 

Axial 

Stress 

Mean  Normal 
Stress 

Principal  Stress 
Difference 

% 

% 

% 

psi 

kPa 

psi 

kPa 

psi 

kPa 

CTCI5_l 

5.0 

-8.8 

13.8 

-12.6 

119.4 

823.1 

50.1 

345.2 

104.0 

716.9 

CTC15_2 

5.2 

-8.5 

13.7 

115.3 

795.5 

on 

337.8 

99.5 

686.6 

CTC15_3 

B 

a 

14.1 

16.5 

113.8 

-11.9 

717.6 

Mean 
(CTC  15) 

5.2 

-8.6 

13.8 

15.9 

118.4 

816.7 

50.1 

345.3 

102.5 

707.0 

0.2 

0.2 

0.2 

0.6 

3.8 

0.5 

18.8 

1.1 

2.6 

17.7 

CTC30_3 

B 

-6.0 

10.5 

31.4 

216.6 

-7.5 

194.3 

1340.1 

85.7 

591.1 

162.9 

1123.6 

CTC30_4 

6.0 

-7.9 

13.9 

31.9 

220.0 

-9.8 

171.1 

1180.2 

78.3 

540.1 

139.2 

960.1 

CTC30_5 

B 

a 

13.1 

31.4 

216.6 

-9.1 

194.2 

1339.2 

85.7 

590.8 

162.8 

1122.6 

Mean 

(CTC30) 

B 

-7.1 

12.5 

31.6 

217.7 

-8.8 

186.5 

1286.5 

83.2 

574.0 

155.0 

1068.8 

0.8 

1.0 

1.8 

0.3 

2.0 

1.2 

13.4 

92.1 

B 

29.4 

13.6 

94.1 

CTC50_lr 

5.1 

-6.2 

11.3 

51.5 

355.2 

-7.2 

274.1 

1890.2 

125.7 

866.9 

222.6 

1535.0 

CTC50__2r 

B 

m 

12.7 

51.3 

353.8 

-9.3 

286.0 

1972.7 

129.5 

893.4 

234.7 

1618.9 

CTC50_3r 

5.6 

-6.6 

12.2 

51.9 

357.9 

-7.6 

268.5 

1851.9 

124.1 

855.9 

216.6 

1494.0 

■ 

B 

IE9 

12.1 

51.6 

355.6 

-8.0 

276.2 

1904.9 

126.4 

872.1 

224.6 

1549.3 

0.2 

0.6 

0.7 

0.3 

2.1 

1.1 

8.9 

61.7 

2.8 

19.3 

B 

63.6 

CTC80_1 

B 

B 

11.4 

76.4 

526.9 

-5.8 

339.9 

2344.3 

164.2 

1132.7 

263.5 

1817.4 

CTC80_2 

5.5 

-5.9 

11.3 

81.6 

562.8 

-6.3 

364.0 

2510.6 

175.7 

1212.1 

282.4 

1947.8 

CTC80_3 

B 

-5.9 

11.6 

80.3 

553.8 

-6.2 

398.7 

2750.0 

186.4 

1285.9 

318.4 

2196.1 

EM 

5.6 

-5.8 

11.4 

79.4 

547.8 

-6.1 

367.6 

2535.0 

175.5 

1210.2 

288.1 

1987.1 

0.1 

0.1 

0.1 

B 

18.7 

0.2 

29.6 

203.9 

11.1 

76.6 

27.9 

192.4 

*  Standard  Deviation 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  80  psi 


Figure  5.11.  Principal  stress  difference  versus  principal  strain  difference  for  conventional 
triaxial  compression  tests  of  granular  limestone 


Crushed  Limestone  Type  610 
Triaxial  Compression 
kpa 


Figure  5.12.  Principal  stress  difference  versus  mean  normal  stress  for  conventional  triaxial 
compression  tests  of  granular  limestone 
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Crushed  Limestone  Type  610 
Triaxial  Compression  at  80  psi 
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Figure  5.13.  Mean  normal  stress  versus  volumetric  strain  for  conventional  triaxial 
compression  tests  of  granular  limestone 

Uniaxial  Strain  Tests 

Uniaxial  strain  tests  were  conducted  until  confining  pressures  reached  a  maximum  of 
100  psi  for  two  tests.  Difficulties  with  membrane  leakage  resulted  in  the  maximum  pressure 
being  reduced  to  80  psi  for  the  two  other  UXE  tests.  Each  UXE  test  was  conducted  by 
applying  an  increment  of  axial  load  until  a  slight  increase  in  specimen  diameter  was  detected. 
Confining  pressure  was  then  applied  until  the  specimen  diameter  returned  to  its  original  value. 
These  processes  were  repeated  throughout  the  test  until  the  desired  maximum  confining 
pressure  was  reached. 

A  summary  of  these  tests  is  shown  in  Table  5.4.  Plots  of  mean  normal  stress  versus 
volumetric  strain  are  shown  in  Figure  5.14.  Figure  5.15  shows  a  composite  plot  of  principal 
stress  difference  versus  principal  strain  difference.  Plots  of  principal  stress  difference  versus 
mean  normal  stress  are  shown  in  Figure  5.16.  The  individual  plots  of  test  data  from  the 
uniaxial  strain  tests  are  shown  in  Appendix  F. 
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Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  5.14.  Mean  normal  stress  versus  principal  stress  difference  for  uniaxial  strain  tests 


Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Principal  Strain  Difference,  % 

Figure  5.15.  Principal  stress  difference  versus  principal  strain  difference  for  uniaxial  strain 
tests 
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Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 
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Figure  5.16.  Principal  stress  difference  versus  mean  normal  stress  for  uniaxial  strain  tests 


Table  5.4.  Summary  of  Peak  Stress  Results  from  Uniaxial  Strain  Tests 


Strain 

Confining 

Pressure 

Volumetric 

Strain 

Axial 

Stress 

Mean  Normal 
Stress 

Principal 

Stress 

Difference 

Axial 

Radial 

Difference 

% 

% 

% 

psi 

kPa 

% 

psi 

kPa 

psi 

kPa 

psi 

kPa 

UXE1 

0.8 

0.1 

0.7 

80.1 

552.4 

1.0 

175.3 

1209.3 

111.9 

771.4 

95.2 

656.8 

UXE2 

1.0 

0.0 

1.0 

80.0 

551.8 

1.0 

171.0 

1179.7 

110.4 

761.1 

91.0 

627.9 

UXE3 

1.0 

0.0 

0.9 

80.0 

551.8 

1.0 

171.5 

1182.6 

110.5 

762.0 

91.5 

630.8 

Mean 

0.9 

0.0 

0.9 

80.0 

552.0 

1.0 

172.6 

1190.5 

110.9 

764.8 

92.6 

638.5 

a* 

0.1 

0.0 

0.1 

0.1 

0.4 

0.0 

m 

16.3 

0.8 

□ 

2.3 

15.9 

*  Standard  Deviation 
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Hydrostatic  Compression  Tests 

Hydrostatic  Compression  Tests  were  conducted  until  confining  pressures  reached  a 
maximum  of  100  psi  (689.5  kPa)  for  four  tests.  The  primary  reason  for  conducting  this  type 
of  test  is  to  define  the  response  of  the  material  in  a  zero  shear  or  pure  normal  stress 
environment.  This  hydrostatic  state  of  stress  produces  strains  that  are  totally  decoupled  from 
any  deviatoric  shear.  The  data  from  this  test  is  used  to  define  the  normal  consolidation  and 
critical  state  parameters  for  the  material.  This  test  is  not  necessarily  representative  of  any 
condition  that  exists.  These  tests  were  conducted  by  applying  an  all  around  pressure  in  the 
CTC  test  device  until  a  desired  maximum  pressure  was  reached.  Difficulties  with  membrane 
leakage  resulted  in  the  maximum  pressure  being  reduced  to  80  psi  (551.6  kPa)  for  one  HC 
test.  A  summary  of  the  peak  stress  results  from  the  HC  tests  is  shown  in  Table  5.5.  Plots  of 
mean  normal  stress  versus  volumetric  strain  are  shown  in  Figure  5.17.  The  individual  plots  of 
test  data  from  the  hydrostatic  compression  tests  are  shown  in  Appendix  D. 


Table  5.5.  Summary  of  Peak  Stress  Results  from  HC  Tests 


Axial 

Strain 

Radial 

Strain 

Hydrostatic 

Pressure 

Volumetric 

Strain 

% 

% 

psi 

kPa 

% 

HC  100 

HC100_1 

0.29 

0.51 

103.9 

716.6 

1.31 

HC100_2 

0.33 

0.49 

99.4 

685.5 

1.30 

HC100_3 

0.36 

0.56 

99.3 

684.9 

1.49 

HC100_4 

0.27 

0.46 

93.4 

644.2 

1.20 

Mean 

0.31 

0.51 

99.0 

682.8 

1.33 

<7* 

0.04 

0.04 

4.3 

29.7 

0.12 

HC  80 

HC80_1 

0.25 

0.50 

79.4 

547.6 

1.15 

*  Standard  Deviation 
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Crushed  Limestone  Type  610 
Hydrostatic  Compression 


Figure  5.17.  Mean  normal  stress  versus  volumetric  strain  for  hydrostatic  compression  tests  of 
granular  limestone 
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DETERMINATION  OF  ABAQUS  DRUCKER-PRAGER 
MODEL  PARAMETERS 

The  values  for  (5  (slope)  and  d  (y  intercept)  for  the  DP  model  were  determined  from  a 
composite  plot  of  the  failure  points  for  the  15,  30,  50,  and  80  psi  (103.4,  206.8,  344.7,  551.6 
kPa)  conventional  triaxial  compression  tests.  The  elliptical  cap  location  is  determined  from 
the  plastic  volume  change  of  a  hydrostatic  compression  test.  The  granular  limestone  material 
tested  was  very  dense  and  strong  when  compared  to  the  types  of  materials  that  the  Drucker- 
Prager  model  was  originally  intended  to  represent.  The  high  density  and  strength  of  the 
material  is  attributed  to  the  high  level  of  compactive  effort  used  in  fabricating  the  specimens 
and  placing  this  material  in  the  field.  As  a  result  of,  the  hydrostatic  stress  regime  under  which 
one  would  see  plastic  volume  change  (i.e.  cap  location)  is  much  higher  than  the  service  loads 
that  even  aircraft  pavements  would  ever  see.  In  essence,  this  reduces  the  cap  model’s 
operation  back  to  a  simpler  two-parameter  friction  model.  Figure  5.18  shows  the  failure 
points  from  the  CTC  test  and  the  hydrostatic  compression  test  plotted  on  a  q  (stress  difference) 
versus  log p  (mean  normal  stress)  space.  Figure  5.19  shows  the  composite  stress  strain 
response  for  the  material.  Figure  6.20  shows  that  the  reference  stress  (virgin  loading)  line  is 
beyond  the  line  bounding  the  maximum  void  ratios  for  the  CTC  test  results.  This  supports  the 
conclusion  that  the  cap  for  this  material  lies  totally  outside  the  stress  regime  of  interest  in  this 
research. 

A  value  of  58.6  0  for  p  and  1 1.25  psi  (77.5  kPa)  for  d  were  calculated  from  the  test 
data.  A  value  of  26,000  psi  (179.3  MPa)  was  calculated  for  the  shear  modulus  of  the  base 
course  as  shown  in  Figure  5.19.  From  the  plot  of  principal  stress  difference  versus  principal 
strain  difference  an  initial  tangent  slope  of  52,000  psi  (358.5  MPa)  was  determined.  This 
value  of  shear  modulus,  G,  is  used  in  calibrating  both  the  Drucker-Prager  and  WES 
Multimechanical  models.  This  initial  shear  modulus  value  is  within  the  normal  valid  range  for 
granular  limestone  base  course  materials  as  seen  in  many  flexible  pavements  (Ulidtz,  1998). 
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Failure  Surface 
Crushed  Limestone 


kpa 
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Mean  Normal  Stress,  psl 
■  Failure  Data  —Failure  Surface 


Figure  5.18.  Failure  surface  for  crushed  limestone  base  course  material 


Triaxial  Tests  of  Granular  Materia! 


Principal  Strain  Difference,  % 

Figure  5.19.  Composite  plot  of  initial  portion  of  principal  stress  difference  versus  principal 
strain  difference  showing  shear  modulus 
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DETERMINATION  OF  WES  MULTIMECHANICAL 
MODEL  PARAMETERS 

The  global  parameters  were  established  using  the  following  methods.  A  summary  of 
the  values  is  shown  in  Table  5.6.  A  stand-alone  version  of  the  WES  MM  model  called 
MVIEWER  was  written  to  aid  in  determining  those  parameters  that  require  trial-and-error 
methods.  MVIEWER  provides  the  analyst  with  a  PC  compatible  platform  to  simulate 
laboratory  tests  relatively  easily.  A  discussion  of  the  MVIEWER  program  and  its  application 
is  presented  in  Appendix  G. 

Strength  Parameters 

The  value  for  friction  angle,  <j),  was  based  on  the  15-psi  (103.4-kPa)  triaxial 
compression  tests  and  the  unconfined  compression  tests.  A  value  of  48  degrees  was 
determined  from  the  tests  using  the  conventional  Mohr’s  Circle  technique  described  in 
Appendix  H.  The  value  for  cohesion  was  selected  based  on  fitting  the  model  to  unconfined 
compression  tests  and  the  15-psi  (103.4-kPa)  CTC  test  data. 

Stiffness  Parameters 

The  bulk  modulus,  K,  is  the  slope  of  the  Mean  Normal  Stress  versus  Volumetric 
Strain  curve  and  was  determined  from  a  hydrostatic  compression  test.  The  shear  modulus  was 
determined  from  a  plot  of  shear  stress  versus  shear  strain.  The  slope  of  the  initial  portions  of 
the  curves  should  be  equal  to  twice  the  shear  modulus,  G.  The  value  picked  for  G  was  more 
than  1.5  times  greater  than  the  K,  thus  the  value  of  Poisson’s  ratio  is  slightly  negative. 
Traditional  engineering  practice  would  consider  reasonable  values  to  range  between  0  and  0.5. 
Theoretically  admissible  values  range  between  -1  and  0.5.  This  material  is  nonlinear, 
anisotropic,  and  plastic  from  early  loading  such  that  the  classical  concept  of  Poisson’s  ratio  is 
not  truly  applicable.  The  value  for  G  represents  the  rate  at  which  shear  stress  accumulates  for 
a  given  amount  a  shear  stress.  A  very  stiff  initial  modulus  is  desirable  so  that  plastic  behavior 
under  relatively  low  stresses  can  be  simulated  using  an  elastic-plastic  model. 
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Other  Parameters 


The  e-  log  p  curve  from  a  hydrostatic  compression  test  should  normally  be  used  to 
obtain  a  slope  for  the  normal  consolidation  line,  NCL.  The  intercept  for  the  line  is  called  the 
hydrostatic  intercept.  The  hydrostatic  compression  tests  did  not  produce  an  overwhelming 
plastic  response,  as  is  the  case  with  clay  or  other  fine-grained  soils.  The  response  had  a 
significant  elastic  component  since  the  preconsolidation  pressure  of  the  base  course  material 
was  not  reached  during  the  test.  Another  approach  was  used  to  produce  an  e-log  p  relation  for 
the  limestone  aggregate. 

The  approach  was  to  plot  the  maximum  void  ratio  achieved  during  the  triaxial 
compression  tests  against  the  logarithm  of  the  mean  stress  associated  with  the  maximum  void 
ratio.  A  plot  of  the  16  data  points  is  shown  in  Figure  5.20.  The  model  line  (gray  lower  line) 
was  drawn  to  provide  an  upper  bound  to  the  data.  The  line  used  for  the  e-log  p  relation  was 
drawn  parallel  to  the  bounding  line  but  with  a  slightly  higher  intercept  of  0.7  psi.  The 
reciprocal  of  Cc  had  a  magnitude  of  8.685  The  values  provided  a  reference  stress,  Pe,  which 
was  used  to  normalize  the  mean  stress  in  the  model. 

The  dilatancy  factor  rate  is  a  scaling  factor  for  shear-volume  coupling.  In  CSSM  only 
one  factor  is  used  to  control  shear-volume  coupling.  The  shear-volume-coupling  factor,  Me, 
is  the  ratio  of  shear  stress  to  mean  stress  under  constant  volume.  Me  should  have  a  value  of 
about  1.8  based  on  the  triaxial  tests  but  this  magnitude  produced  contraction  in  the  model 
during  shearing.  All  of  the  test  specimens  dilated  during  shear  and  the  value  of  Me  was 
adjusted  to  correctly  simulate  the  volume  change  during  shear. 

The  parameters  for  over  consolidation  (OC)  factor  and  phi  ratio  are  intended  to  reduce 
the  strength  of  the  material  as  a  function  of  confining  stress.  The  strength  parameters  for  the 
base  course  aggregate  are  shown  in  Table  5.7.  The  OC  factor  and  phi  ratio  provides  a  function 
to  reduce  phi  as  the  mean  stress  increases.  The  values  adopted  were  chosen  by  trial  and  error. 
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void  ratio 


e  -  log  p 


Figure  5.20.  Plot  used  to  determine  NCL  relationship  for  granular  limestone  material 


Table  5.6.  Global  Properties  for  Granular  Limestone 


PROPERTY 

MAGNITUDE 

BASIS 

Phi 

48  degrees 

tests  at  15  psi  confining  pressure 

Cohesion 

0.25  psi 

Unconfined  Compression 

Bulk  Modulus 

10000  psi 

Hydrostatic  Compression 

Shear  Modulus 

26000  psi 

Plot  of  shear  stress  vs.  shear  strain 

Phi  Ratio 

0.50 

adjust  to  CTC  yield  data 

Hydrostatic 

Intercept  Fh 

0.70  psi 

e  -  log  p  curve 

(or  emax  -  log  p  from  shear  tests) 

Reciprocal  of  Cc 

8.685 

e  -  log  p  curve 

(or  emax  -  log  p  from  shear  tests) 

Shear-Volume 

Factor  MC 

0.72 

adjust  to  volume  change  data 

OC  Factor 

1.80 

adjust  to  yield  data 

Dilatancy  Rate 

Factor 

1.00 

Set  to  unity  as  CSSM  convention 
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Table  5.7.  Strength  Parameters  by  Confining  Stress  for  Granular  Limestone 


Confining  Stress  (psi) 

Mean  Stress  at  Max  Q 

Cohesion  (psi) 

Friction  Angle 

15 

40 

2.1 

48.2 

30 

85 

1.9 

44.5 

50 

124 

1.8 

42.7 

80 

186 

1.6 

39.8 

Mechanism  Parameters 

The  mechanism  parameters  are  shown  in  Table  5.8.  The  underlying  philosophy  for 
obtaining  model  parameters  should  follow  the  approach  outlined  in  the  earlier  section  entitled 
“Calibrating  the  model  -  General  Approach.”  In  fact,  the  calibrations  were  done  in  an 
informal  manner  as  the  more  subtle  features  of  the  model  were  being  discovered  through  the 
act  of  calibration  itself. 


Table  5.8.  Mechanism  Properties  for  Granular  Limestone 


Mechanism 

i 

2 

3 

4 

Phi  Fraction 

0.350 

0.420 

0.820 

0.88 

Mean  Stress  Fraction 

0.900 

0.770 

0.380 

0.48 

Shear  Stiffness  Distribution 

0.702 

0.148 

0.058 

0.0042 

Compression  Limit 

0.018 

0.9 

1.00 

1.00 

Volumetric  Stiffness  Distribution 

0.565 

0.38 

0.02 

0.035 

The  mechanism  parameters  were  adjusted  through  trial  and  error  with  the  MVIEWER 
program  to  conform  to  the  conventional  triaxial  test  data.  The  following  guidelines  proved 
helpful  in  assigning  values  to  the  parameters. 

1 .  Set  the  mechanism  strength  to  yield  for  the  1st  mechanism,  then  2nd,  3rd  and  allow  4th 
mechanism  to  not  yield  at  all.  The  phi  factor  was  used  to  achieve  control  strength. 

2.  Use  the  PFact  to  alter  strength  by  limiting  the  mean  stress  seen  by  each  mechanism. 
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3.  Use  the  shear  ratio  to  adjust  the  stiffness  of  the  shear  mechanism.  The  1st  to  yield  should 
be  stiffest.  The  last  mechanism  to  yield  should  have  the  lowest  stiffness. 

4.  Adjust  Me  to  provide  a  reasonable  volume  change  during  shearing. 

5.  Start  the  calibration  with  bulk  ratio  equal  to  0.25  for  all  mechanisms  and  Hlimit  set  to  unity 
for  all  limits.  Then,  lower  Hlimit  on  1st  mechanism  until  a  effect  is  achieved. 

6.  Special  consideration  is  needed  to  match  the  unconfined  compressive  test  results.  The 
amount  of  mean  stress  applied  to  a  mechanism  should  be  adjusted  with  the  parameter  PFact. 

7.  The  strongest  mechanism  (mechanism  4)  should  be  adjusted  with  PFact  such  that  the 
mechanism  fails  in  unconfined  compression  and  does  not  fail  in  confined  triaxial  compression 
tests. 

APPLICATION  OF  MVIEWER 

The  primary  purpose  of  the  MVIEWER  program  is  to  provide  the  analyst  with  the 
capability  of  easily  evaluating  the  effect  of  changes  in  input  parameter  s  on  the  stress  strain 
response  of  the  model.  The  following  section  demonstrates  this  feature  of  the  MVIEWER  as 
the  program  is  used  to  investigate  the  sensitivity  of  the  stress  strain  response  to  some  of  the 
model  parameters.  As  discussed  earlier,  there  are  ten  global  properties  and  twenty  mechanism 
specific  properties  in  the  WES  MM  model.  The  material  constants  shown  in  Tables  5.6 
through  5.8  contain  the  parameters  used  in  the  FEM  analytical  studies  of  the  laboratory  and 
field  tests. 

Changes  in  the  global  material  properties  have  an  effect  on  the  response  of  all  four 
mechanisms.  The  major  parameters  that  effect  the  shear  strength  are  the  cohesion  (c)  and 
friction  angle  ((f>).  Changes  in  the  shear  modulus  (G)  and  bulk  modulus  (K)  effect  the  stiffness 
of  the  model  response.  The  remaining  global  parameters  are  used  to  adjust  the  model’s 
dilatancy,  and  hydrostatic  response.  When  the  shear  modulus  is  increased,  the  response  of  the 
entire  constitutive  model  stiffens.  In  Figure  5.21,  G  has  been  changed  from  30,000  psi  (206.8 
MPa)  to  60,000  psi  (413.6  MPa). 
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Similar  behavior  occurs  as  G  is  decreased  from  30,000  psi  (206.8  MPa)  to  15,000  psi 
(103.4  MPa).  The  model  response  is  softened  as  shown  in  Figure  5.22.  Changes  in  the  bulk 
modulus,  k,  have  a  smaller  effect  on  the  stress  strain  response  of  the  model  as  shown  in  Figure 
5.23.  As  one  would  expect,  the  K  term  primarily  effects  the  volumetric  strain  response  of  the 
model. 

The  strength  parameters,  C  and  <f>,  have  a  pronounced  effect  on  the  occurrence  of 
yielding  in  the  WES  MM  model.  The  effect  of  increasing  <|)  by  only  10%  (from  48°  to  52.8°  ) 
is  shown  in  Figure  5.24. 

The  granular  limestone  material  used  in  this  study  had  a  very  low  cohesion,  and  large 
changes  in  the  cohesion  parameter  C  had  a  much  less  pronounced  effect  on  yield  than  did 
changes  in  <J>.  The  effect  of  increasing  C  by  10  times  is  picked  up  at  only  the  higher  stress 
level  and  is  shown  in  Figure  5.25. 

The  PHIRATIO  and  DECAY  parameters  are  used  to  adjust  the  friction  angle  as  a 
function  of  mean  stress  and  the  degree  of  dilatancy  experienced.  Figure  5.26  shows  the  small 
effect  that  a  change  in  PHIRATIO  has  on  the  stress  strain  response  of  this  material. 

Figure  5.27  shows  the  effect  of  a  change  in  the  DECAY  parameter  on  the  stress  strain 
response  of  the  WES  MM  model  for  this  type  of  material.  The  effects  of  DECAY  and 
PHIRATIO  are  seen  primarily  in  the  response  of  the  model  at  higher  strain  levels  when 
dilatancy  has  begun  to  occur. 

The  remaining  global  parameters  are  related  to  the  hydrostatic  response  (Me)  of  the 
model  and  dilatancy  scaling  (y).  The  Me  parameter  is  used  to  normalize  the  hydrostatic  stress 
to  the  reference  stress  (Pe)  from  the  Normal  Consolidation  Line.  The  dilatancy  scaling  factor, 
y,  is  left  at  unity  for  the  material  used  in  this  study.  The  twenty  mechanism  specific 
parameters  effect  the  shape  of  the  stress  strain  curve  in  the  same  manner  their  respective 
global  counterparts. 
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Stress  Strain  Curve 


Plot 


Figure  5.22.  Stress  strain  response  with  G=30,0 
psi  (103.4  MPa)  (lower  line) 
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psi  (206.8  MPa)  (upper  line)  and  15,000 


Figure  5.23.  Stress  strain  response  with  K=20,000  psi  (137.9  MPa)  (upper  line),  10,000  psi 
(68.9  MPa)  (middle  line),  and  5,000  psi  (34.3  MPa)  (lower  line) 
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Figure  5.24.  Stress  strain  response  with  <|>=480  (lower  line)  and  <f>=52.6°  (upper  line) 


Figure  5.25.  Stress  strain  response  with  C  =  0.25  (lower  line)  and  C  =  2.5  (upper  line) 
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Figure  5.26.  Stress  strain  response  with  PHIRATIO=0.  5  (lower  line)  and  PHIRATIO=0.75 
(upper  line) 


CHAPTER  6:  MODEL  VERIFICATION 


VERIFICATION  ANALYSES 

In  order  to  verify  that  the  parameters  obtained  in  the  calibration  procedures  provide 
adequate  response  predictions  it  was  necessary  to  conduct  analytical  simulations  of  selected 
laboratory  and  field  tests  using  the  ABAQUS  finite  element  code.  The  predicted  response  was 
compared  to  the  measured  response  obtained  during  the  tests  to  provide  an  indication  of  the 
accuracy  of  the  model  calibrations.  Both  constitutive  models,  Drucker  Prager  and  WES  MM, 
were  used  in  the  laboratory  verification  analyses.  Only  the  WES  model  was  used  in  the  field 
test  analyses. 

ABAQUS  ISSUES 

The  ABAQUS  user  must  enable  some  special  features  and  change  certain  defaults  to 
obtain  a  solution  of  a  non-linear  problem  involving  a  frictional  material.  Three  keywords  are 
used,  *STEP,  *STATIC,  and  ^CONTROL.  The  use  of  the  terms  time  step  and  load  step  are 
used  interchangeably  in  static  analysis  problems  with  ABAQUS.  Time  is  used  as  the  arbitrary 
index  upon  which  loads  are  incremented  to  arrive  at  a  solution  in  all  ABAQUS  runs  in  this 
study.  Each  of  the  examples  comes  from  the  input  file  provided  in  Appendix  A. 

The  *STEP  keyword  should  have  NLGEOM,  EXTRAPOLATION,  UNSYMM,  and 
INC  features  considered.  The  NLGEOM  enables  large  deformation  features.  The 
EXTRAPOLATION  feature  may  be  set  to  EXTRAPOLATION=NO  to  suppress  extrapolation 
of  the  strain  increment  to  the  next  increment.  The  EXTRAPOLATION  =YES  is  the  default. 
The  INC  switch  should  be  used  to  raise  the  number  of  increments  above  the  default  value  of 
10.  Finally,  the  UNSYMM=YES  should  be  enabled  so  that  the  entire  stiffness  matrix  is  used. 
ABAQUS  uses  a  symmetric  matrix  as  default.  Frictional  materials  have  an  unsymmetrical 
stiffness  matrix.  A  symmetric  approximation  may  work,  but  convergence  should  be  made 
easier  if  the  unsymmetrical  matrix  is  used.  Once  NLGEOM  is  turned  on  it  remains  on. 
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♦STEP,  EXTRAPOLATION=NO,  INC=100,  UNSYMM=YES,  NLGEOM 

The  ^STATIC  keyword  has  four  time  parameters  on  one  line:  initial  time  increment, 
total  step  length,  minimum  time  increment  allowable,  maximum  time  increment.  Automatic 
time  incrementation  is  used  whenever  possible.  The  time  has  no  actual  units  of  clock  time 
unless  some  time  dependant  phenomenon  is  being  considered.  ABAQUS  generates  a 
warning  about  this  situation  on  each  static  step.  The  default  for  minimum  time  increment  is 
1.0E-5  time  the  size  of  the  step.  The  static  parameters  must  be  set  each  time  the  *STATIC  is 
used. 

♦STATIC 

0.1,  1.0,  l.E-8,  0.25 

The  *CONTROL  keyword  has  several  important  features.  The  *  CONTROL  keyword 
can  enable  a  line  search  that  is  particularly  important  during  reversals  of  strain  or  stress. 
Enable  the  line  search  and  allow  4  to  6  line  search  iterations.  ABAQUS  uses  a  displacement 
criterion  in  addition  to  a  force  residual  criterion  to  determine  equilibrium.  The  displacement 
criterion  should  be  turned  off  with  the  *CONTROL  keyword.  Finally  the  allowable  number 
of  attempts  to  reach  equilibrium  should  be  increased  above  the  default  values.  The 
♦CONTROL  keyword  can  be  use  to  increase  the  time  incrementation  parameters  to  allow 
more  attempts  to  reach  equilibrium  before  cutting  a  time  increment  and  allow  more  attempts 
to  reach  equilibrium  in  general.  The  defaults  are  too  small  for  non-linear  analyses.  The  control 
options  are  set  once  and  stay  in  force  on  all  subsequent  steps  unless  changed.  When  the 
FIELD  parameter  is  set  to  a  DISPLACEMENT  value  it  affects  the  tolerance  of  residual  force 
that  is  allowed  for  an  increment  to  converge.  The  LINE  SEARCH  parameter  is  used  to  set  the 
solution  technique  to  a  more  robust  form  than  that  used  for  linear  problems.  The  value  is 
basically  a  switch  that  enables  the  line  search  algorithm,  if  the  value  is  anything  other  than 
zero.  The  TIME  INCREMENTATION  parameter  is  set  to  enable  the  automatic  time 
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incrementation  procedures  in  ABAQUS  to  be  changed  to  obtain  solutions  to  highly  non-linear 
problems  (Meade,  1997). 

♦CONTROLS,  PARAMETERS=FIELD,  FIELD=DISPLACEMENT 
0.07,  1.0,  1.0 

(ratio  of  largest  residual  force  to  average  model  force,  displacement  criteria  switches) 
♦CONTROLS,  PARAMETERS=LINE  SEARCH 
6  (enables  line  search  to  be  performed) 

♦CONTROLS,  PARAMETERS=TIME  INCREMENTATION 
12,  18,21,50,  15,,,  15,  ,6 
(incrementation  cut  back  factors) 

These  are  general  findings  and  are  not  absolute  since  each  and  every  new  non-linear 
FEM  model  or  mesh  may  require  different  techniques  to  reach  a  convergent  solution.  Simple 
changes  in  material  properties  may  create  very  complex  model  behaviors  that  become  very 
difficult  to  solve  numerically. 

SIMULATION  OF  LABORATORY  TESTS 
Analytical  simulations  of  all  four  levels  of  conventional  triaxial  compression  tests  [15 
psi  (103.4  kPa),  30  psi  (206.8  kPa),  50  psi  (344.7  kPa),  and  80  psi  (551.6  kPa)],  uniaxial  strain 
tests,  hydrostatic  compression  tests,  and  unconfined  compression  tests  were  conducted  using 
both  models.  The  WES  MM  model  was  used  to  simulate  a  repeated  load  conventional  triaxial 
compression  test  at  low  strain  levels. 

Conventional  Triaxial  Compression  Tests 
Simulations  of  the  4  levels  of  conventional  triaxial  compression  test  were  performed 
using  both  the  DP  model  and  the  WES  MM  model.  The  stress-strain  results  of  the  simulations 
using  both  the  DP  model  and  the  WES  MM  model  are  shown  along  with  the  test  results  in 
Figures  6. 1  through  6.4. 


The  DP  predicted  failure  surface  was  very  close  to  that  obtained  from  the  laboratory 
with  a  predicted  p  of  58.2  °  and  ad  oil  psi  (82.74  kPa)  compared  with  p  =  58.6  °and 
d  =  1 1.25  psi  (130.01  kPa)  from  the  laboratory  data.  The  failure  points  from  the  DP 
simulations  are  shown  in  Figure  6.5.  The  DP  model  behaves  in  a  manner  consistent  with  a 
classical  elastic-plastic  formulation.  The  pre-yield  behavior  is  characterized  by  elastic 
response  followed  by  plastic  response  when  the  yield  stress  is  exceeded.  The  DP  model  was 
calibrated  using  the  procedures  outlined  by  the  ABAQUS  user  documentation  as  described  in 
Chapter  5.  The  model  under  predicts  maximum  stress  in  all  4  simulated  tests.  The  predicted 
stress  strain  behavior  of  the  granular  material  was  quite  different  from  the  response  measured 
in  the  tests.  One  of  the  shortcomings  of  the  DP  model  is  its  inability  to  adequately  capture  the 
response  of  granular  materials  at  low  stress  levels.  Until  the  yield  point  is  reached,  purely 
elastic  recoverable  strain  is  incurred  due  to  load  application.  The  material  will  appear  to  be 
much  stiffer  in  DP  model  predictions  prior  to  yield  than  that  seen  in  tests  since  only  elastic 
behavior  is  modeled  prior  to  the  yield  point. 

The  preconsolidation  pressure  for  the  material  is  well  beyond  the  stress  levels  of 
interest  in  most  pavements.  This  resulted  in  the  cap  portion  of  the  DP  model  not  coming  into 
play  for  any  of  the  simulations,  which  produced  non-dilative  response  predictions.  The  non- 
dilative  behavior  of  the  DP  model  can  be  see  in  Figure  6.6,  where  the  mean  normal  stress  is 
plotted  versus  volumetric  strain. 

The  agreement  of  the  WES  MM  model  predictions  with  the  test  data  is  best  for  the  30- 
psi  (206.8  kPa)  test.  The  stress-strain  response  for  all  four  CTC  tests  are  very  good.  The  basic 
shape  of  the  curve  and  the  maximum  stress  level  reached  is  very  close  for  all  except  the  80-psi 
(551.6  kPa)  test.  The  accuracy  of  the  80-psi  (551.6  kPa)  test  prediction  was  sacrificed  to 
achieve  a  closer  fit  at  lower  stress  levels  (i.e.,  stress  levels  closer  to  those  expected  in  field 
tests  to  be  discussed  later  in  this  chapter).  The  predicted  failure  surface  of  the  WES  MM 
model,  in  a  principal  stress  difference  versus  mean  normal  stress  space,  was  very  close  to  that 
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obtained  from  the  laboratory  with  a  predicted  friction  angle,  p,  of  60  °  and  a  cohesion,  d,  of 
12  psi  (82.74  kPa)  compared  with  P  =  58.6  °  and  d=  11.25  psi  (130.01  kPa)  from  the 
laboratory  test  data.  A  plot  of  the  WES  MM  failure  data  is  shown  in  Figure  6.7.  This  plot  is 
very  close  to  the  failure  surface  plotted  from  the  laboratory  test  data.  The  ability  of  the  WES 
MM  model  to  change  friction  angle  with  increasing  mean  normal  stress  is  primarily 
responsible  for  the  small  differences. 

The  post  yield  shear  dilatant  behavior  of  the  WES  MM  model  is  also  demonstrated  in 
the  composite  plot  of  mean  normal  stress  versus  volumetric  strain  shown  in  Figure  6.8.  The 
breakpoints  (changes  in  slope  of  the  stress-strain  curve)  in  the  WES  MM  model  are  also  very 
evident  in  the  response  shown  in  Figure  6.8. 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  15  psi  (103.4  kPa) 
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Figure  6. 1 .  Composite  plot  of  principal  stress  difference  versus  principal  strain  difference  for 
15  psi  (103.4  kPa)  test 
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Crushed  Limestone  Type  610 
Triaxial  Compression  at  30  psi  (206.8  kPa) 


Figure  6.2.  Composite  plot  of  principal  stress  difference  versus  principal  strain  difference  for 
30  psi  (206.8  kPa)  test 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  50  psi  (344.7  kPa) 


Figure  6.3.  Composite  plot  of  principal  stress  difference  versus  principal  strain  difference  for 
50  psi  (344.7  kPa)  test 


108 


Crushed  Limestone  Type  610 
Trlaxial  Compression  at  80  psi  (551.6  kPa) 


Figure  6.4.  Composite  plot  of  principal  stress  difference  versus  principal  strain  difference  for 
80  psi  (551.6  kPa)  test 
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Figure  6.5.  Predicted  failure  surface  for  Drucker-Prager  model  compared  with  test  results 
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Crushed  Limestone  Type  610 
Drucker  Prager  Model 


Figure  6.6.  Composite  plot  of  mean  normal  stress  versus  volumetric  strain  for  DP  predictions 
of  CTC  tests 
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Figure  6.7.  Predicted  failure  surface  for  WES  MM  model  compared  with  test  results 
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Crushed  Limestone  Type  610 
WES  MM  Model 


Figure  6.8.  Composite  plot  of  mean  normal  stress  versus  volumetric  strain  for  WES  MM 
predictions  of  CTC  tests 


Uniaxial  Strain  Tests 

Simulations  of  the  uniaxial  strain  tests  (UXE)  were  performed  using  both  the  DP 
model  and  the  WES  MM  model.  The  stress-strain  results  of  the  UXE  simulations  using  the 
DP  model  and  the  WES  MM  model  are  shown  with  the  test  results  in  Figures  6.9  through 
6.11.  The  WES  MM  model  demonstrates  the  ability  to  predict  the  stress  path  required  to 
maintain  uniaxial  strain  conditions  through  loading  and  unloading.  The  existence  of  a  residual 
(locked  in)  stress  is  typical  for  this  kind  of  test,  and  can  be  seen  in  the  WES  MM  model  as 
well  as  the  test  data  shown  in  Figure  6.9.  The  DP  model  is  unable  to  capture  the  shearing 
response  that  creates  the  “locked  in”  stress  after  unloading.  This  can  be  attributed  to  the  fact 
that  the  UXE  stress  path  never  intersects  the  DP  failure  surface.  That  leaves  only  elastic 
response  for  the  DP  model  during  this  test.  The  WES  MM  model  is  able  to  capture  this 
behavior  through  the  separation  of  hydrostatic  and  shear  response  in  the  HYDROS  and 
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Ammos  routines  described  in  Chapter  4.  The  plot  of  mean  normal  stress  versus  volumetric 
strain  shown  in  Figure  6. 1 1  demonstrates  the  ability  of  the  WES  model  to  capture  the  overall 
stiffness  of  the  material  in  hydrostatic  conditions,  but  its  shortcomings  in  modeling  the  true 
volumetric  stress-strain  response  are  also  evident.  Additional  effort  in  the  calibration  of  the 
hydrostatic  mechanism  parameters  of  the  WES  MM  model  may  well  provide  the  accuracy  of 
response  predictions  missing  in  these  analyses.  Although  this  type  of  test  is  useful  in 
exercising  the  model,  it  is  not  particularly  representative  of  any  real  condition  that  exists  in  a 
loaded  pavement. 


Ousted  Limestone  Type  610 
Ukiaxiai  Strain  Test 

kpa 

0  100  200  300  400  500  600  TOO  800 


Figure  6.9.  Composite  plot  of  principal  stress  difference  versus  mean  normal  stress  for 
uniaxial  strain  tests 
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Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 
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Figure  6.10.  Composite  plot  of  principal  stress  difference  versus  principal  strain  difference  for 
uniaxial  strain  tests 


Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 
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Figure  6. 1 1.  Composite  plot  of  mean  normal  stress  versus  volumetric  strain  for  uniaxial  strain 
tests 
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Hydrostatic  Compression  Tests 


Simulations  of  the  hydrostatic  compression  tests  (HC)  were  performed  using  both  the 
DP  model  and  the  WES  MM  model.  The  stress-strain  results  of  the  HC  simulations  using  the 
DP  model  and  WES  MM  model  are  shown  with  the  test  results  in  Figure  6.12.  The  DP  model 
can  only  produce  a  linear  response  to  hydrostatic  state  of  stress  that  is  a  function  of  the  bulk 
modulus  that  is  fixed  by  selection  of  Young’s  modulus  and  Poisson’s  ratio.  The  response  of 
the  DP  model  is  stiffer  than  the  test  data  under  hydrostatic  conditions.  However,  matching  the 
response  of  the  DP  model  under  CTC  conditions  was  more  crucial  than  matching  hydrostatic 
response.  Therefore  the  Young’s  modulus  and  Poisson’s  ratio  were  selected  with  the  CTC 
tests  as  the  benchmark  test  results.  Such  a  trade  off  in  performance  is  a  shortfall  of  a 
simplistic  model  like  DP.  In  essence,  one  only  has  four  parameters  to  work  with  to  produce 
stiffness  and  yield  that  can  fit  only  a  limited  number  of  situations.  The  WES  MM  model  did 
not  perform  as  well  for  hydrostatic  test  conditions  as  it  did  in  simulations  of  tests  with  a  lot  of 
shear  stress.  The  response  of  the  model  is  almost  purely  linear  and  does  not  exhibit  yield  in 
hydrostatic  test  conditions.  As  was  the  case  with  the  DP  model,  matching  the  response  of  the 
WES  MM  model  under  CTC  conditions  was  more  crucial  than  matching  the  response  of 
hydrostatic  test  conditions. 

Unconfined  Compression  Tests 

Analytical  simulations  of  the  unconfined  compression  tests  (UCC)  were  performed 
using  both  the  DP  model  and  the  WES  MM  model.  The  stress-strain  results  of  the  HC 
simulations  using  the  DP  model  and  WES  MM  model  are  shown  with  the  test  results  in 
Figures  6.13  and  6.14.  The  DP  model  provides  an  acceptable  prediction  of  the  yield  stress, 
however  it  does  only  a  minimal  job  of  modeling  the  overall  stress  strain  response  of  the 
unconfined  tests.  Again,  a  model  like  the  DP  is  just  to  simple  in  nature  to  capture  the  complex 
response  of  an  unconfined  test  of  granular  material.  The  WES  MM  model  did  a  good  job  of 
modeling  the  overall  stress  strain  response  of  the  unconfined  compression  test.  The  WES 


114 


MM  model  also  has  the  ability  to  capture  both  the  dilative  and  softening  response  of  the 
granular  material  under  unconfined  compression.  These  test  also  serve  primarily  as 
calibration  tests  and  an  index  of  a  models  applicability  to  a  wide  range  of  conditions  and  are 
not  representative  of  conditions  that  would  exist  in  the  granular  base  course  of  a  pavement. 

Cyclic  Triaxial  Compression  Tests 

Simulations  of  the  cyclic  triaxial  compression  tests  (CTCR)  were  performed  using  the 
WES  MM  model.  The  analytical  predictions  made  using  the  standard  calibration  presented  in 
Chapter  5  produced  cyclic  behavior  that  was  somewhat  different  form  the  laboratory  tests. 

The  model  did  produce  hysteresis,  but  the  shape  and  size  of  the  hysteresis  loops  at  low  strain 
levels  and  the  magnitude  of  the  strain  at  which  the  cyclic  behavior  began  was  different  from 
the  test  data  as  shown  in  Figure  6.15.  A  modified  calibration  was  completed  (Shown  in  Table 
6. 1  and  Table  6.2),  and  the  cyclic  test  was  rerun  with  the  new  calibration.  The  ability  of  the 
WES  MM  model  to  closely  capture  cyclic  response  was  clearly  demonstrated  in  this  analysis. 
However,  obtaining  this  calibration  required  a  lot  of  iterations  and  intimate  knowledge  of 
model  behavior.  One  of  the  original  goals  of  this  model  development  was  to  produce  a 
constitutive  model  that  would  be  relatively  easy  to  calibrate  from  standard  geotechnical 
laboratory  tests.  This  modified  calibration  also  proved  to  create  numerical  convergence 
problems  with  ABAQUS  when  applied  to  the  field  test  section  FEM  analyses.  When  using  a 
commercial  finite  element  code  like  ABAQUS,  one  does  not  have  access  to  the  source  code 
for  the  finite  element  program.  As  a  result,  when  problems  with  convergence  are  encountered 
and  can  not  be  solved  through  the  use  of  *CONTROL  options,  other  avenues  of  completing  an 
analysis,  such  as  equivalent  alternate  material  model  calibrations,  must  be  considered.  Even 
though  the  modified  cyclic  calibration  produces  excellent  stress-strain  agreement  with  the  test 
data  at  low  stress  levels,  the  amount  of  permanent  strain  accumulated  from  each  cycle  is  very 
close  for  both  calibration.  Given  these  considerations,  the  original  standard  calibration  was 
used  for  all  test  section  analyses. 
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Table  6.1.  Global  Properties  for  Modified  Calibration 


PROPERTY 

MAGNITUDE 

Phi 

48  degrees 

Cohesion 

0.25  psi  (1.72  kPa)  j 

Bulk  Modulus 

100000  psi  (689.5  MPa) 

Shear  Modulus 

200000  psi  (689.5  MPa) 

Phi  Ratio 

0.50 

Hydrostatic  Intercept 

0.70  psi  (4.82  kPa) 

Reciprocal  of  Cc 

8.685 

Shear-volume  Factor  Me 

0.72 

OC  factor 

1.80 

Dilatancy  Rate  Factor 

1.00 

Table  6.2.  Mechanism  Properties  for  Modified  Calibration 


Mechanism 

1 

2 

3 

4 

Phi  Fraction 

0.1 

0.25 

0.6 

0.9 

Mean  Stress  Fraction 

2.2 

0.86 

0.3 

0.35 

Shear  Stiffness  Distribution 

0.49 

0.26 

0.068 

0.011 

Compression  Limit 

0.018 

0.9 

1.00 

1.00 

Volumetric  Stiffness  Distribution 

0.565 

0.38 

0.02 

0.035 
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Crushed  Limestone  Type  610 
Hydrostatic  Compression  to  100  psi 
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Figure  6.12.  Composite  plot  of  mean  normal  stress  versus  volumetric  strain  for  hydrostatic 
compression  tests 
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Figure  6.13.  Composite  plot  of  axial  stress  versus  axial  strain  for  unconfined  compression 
tests 
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Crushed  Limestone  Type  610 
Unconfined  Compression 
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Figure  6. 14.  Composite  plot  of  mean  normal  stress  versus  volumetric  strain  for  unconfmed 
compression  tests 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  50  psi  (344.7  kPa) 
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Figure  6.15.  Comparison  of  FEM  prediction  of  cyclic  response  with  test  data 
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FIELD  TEST  SECTIONS 


Cyclic  tire  loads  of  two  test  items  from  the  selected  test  section  (Webster,  1993)  were 
simulated  with  the  WES  MM  model  being  used  to  define  the  properties  of  the  granular  base 
course  material.  The  remaining  layers  were  modeled  using  linear  elastic  properties. 

General  Test  Section  Description 

Test  sections  were  located  at  the  U.S.  Army  Engineer  Waterways  Experiment  Station 
(WES)  in  Vicksburg,  MS.  They  were  constructed  within  an  aircraft  hangar  so  that  they  could 
be  sheltered  from  rain  and  sun.  The  existing  soil  floor  was  excavated  to  a  depth  of  40  in. 
(1016  mm)  and  the  lean  clay  at  the  bottom  of  the  trench  was  compacted  to  a  CBR  strength 
greater  than  10.  The  bottom  and  sides  of  the  trench  were  lined  with  sheets  of  polyethylene  to 
minimize  drying  of  the  heavy  clay  subgrade  during  traffic  tests. 

The  subgrade  under  all  test  items  consisted  of  heavy  clay  (CH)  material,  according  to 
the  Unified  Soil  Classification  System.  This  material  had  a  liquid  limit  (LL)  of  67  and  a 
plasticity  index  (PI)  of  45.  When  compacted  in  accordance  with  ASTM  D  698  (standard 
Proctor),  it  had  an  optimum  moisture  content  of  23  percent  (by  mass  of  dry  material), 
corresponding  to  a  maximum  dry  density  of  92  lb/ft3  (1475  kg/m3 ).  Compaction  was 
accomplished  in  6-in.  (152-mm)  lifts  with  a  rubber-tired  roller.  The  final  subgrade  surface 
was  smoothed  with  a  vibratory  steel  drum  roller. 

The  base  course  material  was  an  MDOT  type  610,  crushed  limestone.  When 
compacted  in  accordance  with  ASTM  D  1557  (modified  Proctor),  it  had  an  optimum  moisture 
content  of  4.5  percent  (by  mass  of  dry  material),  corresponding  to  a  maximum  dry  density  of 
144  lb/ft3  (2307  kg/m3).  The  base  course  material  was  back-dumped,  spread  with  a  bulldozer 
and  compacted  in  6-in.  (152-mm)  lifts  with  a  vibratory  steel  drum  roller.  The  top  lift  of  the 
base  course  was  also  compacted  with  a  solid  rubber-tired  roller.  All  test  items  were  surfaced 
with  2  in.  (5 1  mm)  of  asphalt  concrete  (50-blow  Marshall  specification).  The  maximum 
aggregate  size  for  the  asphalt  concrete  was  0.5  in  (12.5  mm)  and  the  minimum  Marshall 
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stability  was  1500  lbs.  (6.67  kN).  The  asphalt  surfacing  covered  an  area  50  ft  (15.2  m)  by  224 
ft  (68.3  m).  This  area  included  all  test  items  and  extended  40  ft  (2.2  m)  past  each  end  of  the 
test  section. 

As-Constructed  Properties 

Cross-section  level  readings  taken  during  construction  indicated  that  the  base  layer 
thicknesses  of  all  test  items  were  constructed  to  within  1  in.  (25  mm)  of  design  thickness. 
Thickness  estimates  for  the  asphalt  concrete  were  obtained  from  twelve  core  samples. 
Thicknesses  ranged  from  2.0  in.  (51  mm)  to  2.6  in  (66  mm).  Table  6.3  shows  a  summary  of 
as-constructed  data  for  base  course  and  subgrade. 


Table  6.3.  As-constructed  Properties  for  Subgrade  and  Base  Course. 


Test 

Lane 

Material 

Depth  of 
Measurement 
in.  (mm) 

CBR 

Moisture 

Content 

(%) 

Dry  Density, 
lb/fit3  (kg/m3 ) 

Percent  of 
Maximum 
Density 

1 

base  course 

2  (51) 

1 

1 

1 

1 

1 

1 

4.1 

137.2  (2198) 

96 

i 

subgrade 

10  (254) 

7.1 

26.3 

92.8  (1487) 

102 

i 

16  (406) 

6.9 

26.2 

92.3  (1479) 

101 

i 

22  (559) 

n 

25.9 

93.5  (1498) 

103 

2 

base  course 

2  (51) 

■ 

t 

i 

i 

• 

i 

mm 

136.9  (2194) 

95 

2 

subgrade 

16  (406) 

2.5 

31.4 

86.6(1388) 

98 

2 

22  (559) 

2.7 

30.5 

86.9(1393) 

98 

2 

28  (711) 

2.3 

31.9 

86.0(1378) 

98 

Instrumentation 

Multi-depth  deflectometers  (MDD)  were  used  to  measure  both  recoverable  and 
permanent  deformations.  A  single  MDD  was  installed  in  each  of  test  items  1  and  2  in  both 
traffic  lanes  1  and  2.  These  MDDs  consisted  of  a  support  shaft  and  up  to  four  modules,  each 
of  which  measured  vertical  deflection  at  a  different  depth.  In  plan  view,  the  MDDs  were 
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centered  in  the  test  items.  The  MDDs  were  retrofitted  into  the  pavements  by  using  an  auger  to 
dig  a  vertical  hole,  1.5  in.  (37.5  mm)  in  diameter,  through  the  pavement  system.  The  MDD 
shafts  were  anchored  8  ft  (2.4  m)  below  the  surface  in  order  to  provide  a  motionless  reference 
for  deflection  measurements.  Each  MDD  included  a  module  at  a  depth  just  beneath  the 
asphalt  concrete  layer  and  just  beneath  the  base  course  layer.  Figure  6.16  shows  a  general 
schematic  of  an  MDD. 
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Figure  6.16.  Typical  cross-section  of  MDD  after  installation 
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Traffic 


Test  traffic  was  applied  using  a  single-wheel  test  cart  with  a  wheel  load  of  30,000  lb 
(133  kN)  as  shown  in  Figure  6.17.  The  loaded  tire,  which  was  designed  for  a  C-130  aircraft, 
was  inflated  to  provide  a  contact  pressure  of  68  psi  (468  kPa).  The  contact  area  was  442  in2 
(0.29  m2 ).  The  cart  was  powered  by  the  front  half  of  a  four-wheel  drive  truck  and  was 
equipped  with  an  outrigger  wheel  to  prevent  overturning.  The  load  was  produced  using  lead 
blocks  located  at  the  rear  of  the  cart.  The  test  traffic  was  applied  to  the  pavement  using  the 
powered  cart  in  both  directions  with  the  drivers  manually  operating  and  aligning  the  vehicle  to 
insure  that  proper  load  distributions  were  maintained. 


Figure  6.17.  Loaded  single  wheel  test  cart  with  C- 1 30  tire 
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ABAQUS  FEM  ANALYSIS  OF  TEST  SECTION 

The  ABAQUS  FEM  code  was  used  to  analyze  the  response  of  the  two  test  sections  of 
interest.  All  ABAQUS  computations  were  conducted  on  SGI  ORIGIN  2000  supercomputers. 
Finite  element  model  development  for  ABAQUS  was  accomplished  interactively  on  engineer¬ 
ing  workstations  using  The  MacNeal-Schwendler  Corporation’s  PATRAN  software  incorpo¬ 
rating  an  ABAQUS  application  interface.  PATRAN  was  also  utilized  to  post-process  many  of 
the  results  from  ABAQUS.  A  2-D  static  axisymmetric  analysis  was  performed  using  the  WES 
MM  model  for  the  base  course  and  linear  elastic  properties  for  the  asphalt  and  subgrade 
layers.  The  purpose  for  this  analysis  was  to  demonstrate  the  ability  to  predict  permanent 
deformation  in  a  granular  pavement  layer  using  a  non-linear  elastic-plastic  model. 

Material  Properties 

The  information  available  from  the  test  section  data  did  not  allow  for  direct  calibration 
of  the  asphalt  and  subgrade  properties.  However,  there  was  enough  information  to  arrive  at 
reasonable  values  for  the  elastic  constants:  Young’s  modulus  (E)  and  Poisson’s  ratio  (v). 
Typical  values  for  these  material  constants  can  be  found  in  a  number  of  sources  (Kulhawy, 
Mayne,  1990)  (Ulidtz,  1998)  (Tseng,  1988).  The  MDD  deformation  values  were  used  to  fine- 
tune  these  constants  to  produce  reasonable  values  of  deformation  under  load.  In  essence  the 
MDD  reading  under  load  enabled  a  crude  backcalculation  of  elastic  constants  to  be  performed, 
thus  enabling  the  base  course  layer  to  see  a  stress  state  very  similar  to  the  true  state  of  stress 
under  load.  Table  6.4  gives  the  values  used  for  the  elastic  constants  in  both  test  sections.  The 
values  used  for  these  elastic  constants  in  the  subgrade  material  represent  more  than  just  a 
material  property.  They  are  an  effective  foundation  stiffiiess  that  allows  for  the  subgrade 
material,  lower  supporting  layers  and  far  boundaries  to  be  included  in  a  simplified  system 
model.  The  crushed  limestone  aggregate  base  course  was  modeled  with  the  WES  MM  model 
as  calibrated  in  Chapter  5. 
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Table  6.4.  Material  Properties  Used  for  Asphalt  and  Subgrade  Layers 


Section  ID 

Asphalt  Layer 

Subgrade  Layer 

E,  psi  (MPa) 

V 

E,  psi  (MPa) 

Lane  1-1 

500,000  (3447.5) 

0.35 

18,000(124.1) 

0.35 

Lane  2-1 

500,000  (3447.5) 

0.35 

15,000  (103.4) 

0.35 

FEM  Mesh 

The  two  test  sections  were  modeled  using  standard  4  node  quadrilateral  axisymmetric 
elements  from  the  ABAQUS  element  library  as  shown  in  Figure  6.18.  In  order  to  correctly 
model  the  C-130  wheel  load,  the  load  was  simulated  with  a  surface  pressure  of  68  psi  (468.8 
kPa)  applied  over  a  circular  area  of 442  square  inches  (0.29  m2)  to  produce  a  total  load  of 
30,000  lbs.  (133  kN).  The  nominal  layer  thickness  values  are  shown  in  Table  6.5.  The  total 
depth  of  the  subgrade  was  240  in.  (6096  mm)  yielding  a  total  model  depth  of  not  less  than  20 
feet  (6.25  m).  This  depth  was  similar  to  that  arrived  at  for  similar  analyses  in  the  literature. 
(Bryant  1998)  (Yeh,  1989)  (Barksdale,  1973) 


Table  6.5.  Layer  Thickness  Values 


Section  ID 

Asphalt  Layer 

Base  Course 

Subgrade 

Lane  1-1 

2  in  ( 50.8  mm) 

10  in  (254.0  mm) 

240  in  (6096  mm) 

Lane  2-1 

2  in  ( 50.8  mm) 

18  in  (457.2  mm) 

240  in  (6096  mm) 

The  coordinate  system  used  for  the  analysis  denotes  Y  as  the  vertical  direction  and  X 
as  the  horizontal  direction.  The  meshes  for  both  test  sections  were  fixed  in  the  X  along  the  left 
side  (line  of  symmetry)  and  the  right  hand  side.  The  meshes  were  fixed  against  vertical 
translation  along  the  bottom.  Figures  6.19  a  and  b  show  the  FEM  mesh  for  the  test  sections. 
The  elements  ranged  in  size  from  the  3  in.  by  2  in.  elements  directly  under  the  loaded  area  to 
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the  12  in.  square  elements  in  the  region  farthest  from  the  loaded  area.  The  elements  at  the 
upper  right  hand  comer  and  lower  left  hand  comer  of  the  model  had  aspect  ratios  on  the  order 
of  6: 1  resulting  in  a  few  slender  elements.  Although  these  aspect  ratios  are  fairly  high, 
ABAQUS  will  support  elements  in  this  aspect  ratio  range,  and  the  area  undergoing  load 
application  was  also  very  far  from  the  slender  elements. 

The  objective  of  this  study  was  to  develop  a  response  model  for  granular  pavement 
layers.  Although  that  primarily  involved  the  development  and  implementation  of  a 
constitutive  model,  it  was  necessary  to  apply  said  model  to  the  analysis  of  a  pavement  system. 
For  validation  As  the  analysis  phase  of  the  study  progressed  it  became  obvious  that  obtaining 
a  convergent  solution  with  ABAQUS  was  a  difficult  task  that  changed  with  the  inclusion  of 
small  differences  in  mesh  or  material  model  parameters.  The  FEM  grids  shown  in  Figures 
6. 19  and  6.20  were  arrived  at  through  trial  and  error  attempts  at  defining  the  finest  mesh  that 
would  provide  reasonable  response  while  still  being  able  to  converge  to  a  solution  during  load 
application.  The  smallest  element  in  these  meshes  is  2  in.  by  3  in.  In  the  earliest  meshes 
developed  during  this  study  the  elements  were  1  in.  by  1  in.  under  the  loaded  area.  The  nature 
of  non-linear  plastic  analysis  proved  to  counter  intuitive  to  traditional  finite  element  mesh 
concepts.  In  most  cases  the  finer  a  mesh  is  made  the  easier  and  more  accurate  the  solution 
will  become  until  the  accuracy  reaches  some  asymptotic  value.  That  concept  relies  on  the  fact 
that  a  solution  can  be  obtained.  In  the  case  of  non-linear  elastic-plastic  analysis  there  is  a  limit 
on  the  minimum  element  size  that  can  provide  a  practical  solution  with  reasonable  load  step 
sizes.  As  the  size  of  the  elements  undergoing  plastic  deformation  gets  smaller  it  is  also 
necessary  to  reduce  the  strain  increment  applied  to  those  elements.  If  the  strain  increment  is 
large  relative  to  the  element  size  then  the  entire  element  and  its  neighbors  may  yield  in  one 
increment.  If  that  happens  then  the  solution  procedures  used  in  codes  like  ABAQUS  will  not 
be  able  to  converge. 
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A  convergent  solution  was  not  possible  until  elements  with  minimum  dimensions  of  at 
least  2  inches  (5.08  mm)  were  used  under  the  loaded  area.  Even  then,  extreme  care  was 
required  in  selecting  the  load  step  (strain  increment)  parameters  required  to  obtain  a 
convergent  solution..  The  formulation  of  the  WES  MM  model  produces  breakpoints  in  the 
stress  strain  curve  each  time  one  of  the  four  mechanisms  in  the  model  yields.  In  a  pavement 
section  FEM  model  any  element  in  the  base  course  was  subject  to  yielding  at  up  to  four  levels 
under  load  application  due  to  these  breakpoints.  As  the  load  was  increased  in  each  of  the 
cycles  of  applied  wheel  load,  the  area  of  the  model  under  the  load  would  experience  yield. 

The  difficulty  in  reaching  a  convergent  solution  would  increase  each  time  the  zone  of  plastic 
deformation  would  move  far  enough  through  the  base  course  to  encompass  another  element. 

In  effect,  plastic  behavior  in  the  material  translates  directly  into  increased  difficulty  in 
obtaining  a  convergent  solution.  This  difficulty  is  coupled  with  the  mesh  fineness  to  produce 
a  very  complex  challenge  in  conducting  an  analysis.  The  automatic  time  stepping  methods 
used  in  ABAQUS  enabled  the  user  to  specify  parameters  for  controlling  the  size  of  the  load 
step  (strain  increment).  ABAQUS  could  use  very  small  steps  when  solution  convergence  was 
difficult  to  obtain,  and  then  use  relatively  large  steps  when  the  model  was  more  numerically 
stable.  The  exact  mesh  dimensions  and  load  step  definitions  can  not  be  obtained  through  a 
direct  method,  but  they  are  arrived  at  obtained  through  trial  and  error.  Unlike  many  trial  and 
error  methods  this  procedure  is  basically  a  “GO”  or  “NO  GO”  proposition.  One  must  also 
weight  the  advantages  of  a  very  accurate  material  calibration  at  low  stress  levels  as  compared 
with  the  increased  difficulty  in  obtaining  a  solution.  The  cyclic  calibration  produced  great 
response  predictions  for  calibration  verifications,  but  could  not  be  made  to  converge  with  a 
pavement  test  section  grid.  The  grids  shown  in  Figures  6. 19  and  6.20  did  converge  with  the 
standard  material  calibration  arrived  at  in  Chapter  5  and  the  ABAQUS  model  definitions 
shown  in  Appendix  B. 
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Figure  6.18.  Typical  axisymmetric  FEM  model  of  a  pavement 
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Figure  6.20.  Finite  element  mesh  for  Lane  2-1 
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Results  of  FEM  Analyses  of  Test  Sections 

The  test  sections  were  subjected  to  5  cycles  of  a  simulated  single  C-130  tire  load.  The 
base  course  layer  was  modeled  with  the  WES  MM  model,  while  the  remaining  layers  were 
modeled  as  a  linear  elastic  material.  The  results  of  these  analyses  were  compared  to  MDD 
measurements  from  the  field  test  sections  to  provide  model  validation  and  assessment. 

Figure  6.21  shows  the  deformed  shape  of  Lane  1-1  under  the  5th  load  application.  The 
value  of  deformation  at  the  top  of  the  base  course  was  173  mils  (4.39  mm).  Webster,  1993, 
reported  the  value  of  deformation  at  the  top  of  the  base  course  in  Lane  1-1  to  be  165  mils 
(4. 19  mm)  under  the  5th  load  application.  The  predicted  deformation  under  load  at  the  top  of 
the  subgrade  was  1 13  mils  (2.87)  as  compared  with  a  field  value  of  125  mils  (3.17  mm).  The 
agreement  between  these  values  verifies  the  relative  accuracy  of  the  overall  system  calibration 
for  Lane  1-1. 

In  order  to  validate  the  ability  of  the  WES  MM  model  to  predict  plastic  accumulated 
strain  under  repeated  loads,  the  value  of  permanent  deformation  after  removal  of  the  load  at 
the  top  of  the  base  course,  33  mils  (0.83  mm),  was  determined  from  the  analysis  and  compared 
with  the  field  value  of  40  mils  (1.02  mm)  for  Lane  1-1.  Figure  6.22  shows  the  deformed  shape 
of  Lane  1-1  after  removal  of  the  5th  load. 

Similar  finite  element  calculations  were  made  for  Lane  2-1.  Figure  6.23  shows  the 
deformed  shape  of  Lane  2-1  under  the  5th  load  application.  Figure  6.23  shows  the  deformed 
shape  of  Lane  2-1  after  removal  of  the  5th  load.  The  deformation  at  the  top  of  the  base  course 
was  predicted  to  be  243mils  (6.17  mm).  Webster,  1993,  reported  the  value  of  deformation  at 
the  top  of  the  base  course  in  Lane  2-1  to  be  190  mils  (4.83  mm)  under  the  5th  load  application. 
The  predicted  permanent  deformation  at  the  top  of  the  base  course  was  93  mils  (2.36  mm)  as 
compared  with  the  field  value  of  50  mils  (1 .27  mm).  The  predicted  deformation  under  load  at 
the  top  of  the  subgrade  was  1 16  mils  (2.94  mm)  with  a  field  value  of  145  mils  (3.68  mm). 
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Figures  6.25  and  6.26  show  the  deformation  under  load  at  the  top  of  the  base  course 
and  the  top  of  the  subgrade  from  the  FEM  predictions  and  the  field  measurements  for  Lanes  1- 
1  and  2-1  respectively.  Figures  6.25  and  6.26  also  show  the  FEM  prediction  of  permanent 
deformation  as  a  function  of  load  cycles.  The  agreement  between  the  FEM  predictions  and 
the  field  measurements  are  closer  for  Lane  1-1  than  for  Lane  2-1 .  The  effects  of  the  lower 
strength  subgrade  (3  CBR)  in  Lane  2-1  are  much  more  difficult  to  model  with  linear  elasticity 
than  the  higher  strength  (8  CBR)  subgrade  in  Lane  1-1.  The  lower  strength  subgrades  are  also 
much  harder  to  construct  and  much  more  susceptible  to  environmental  changes,  which 
resulted  in  higher  variability  of  strength  and  stiffness  within  the  test  section  (Webster,  1993). 
The  magnitude  of  these  deformations  is  very  small  when  compared  with  the  overall  cross- 
sectional  dimensions  modeled.  For  instance,  the  deformation  at  the  top  of  the  base  course 
under  load  in  Lane  2-1  was  predicted  to  be  243  mils  (6. 17  mm).  This  was  53  mils  greater  than 
the  field  measurement  and  only  0.3%  of  the  layer  thickness. 


Figure  6.21.  Deformed  shape  (100  X)  under  5th  load  application  for  Lane  1-1 
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Figure  6.22.  Deformed  shape  (100  X)  after  5th  load  application  for  Lane  1-1 


Figure  6.23.  Deformed  shape  (100  X)  under  5th  load  application  for  Lane  2  -1 
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Figure  6.24.  Deformed  shape  (100  X)  after  5th  load  application  for  Lane  2-1 


32 


Vertical  Deformation,  inches 


Lane  1-1 


Figure  6.25.  Vertical  deformation  versus  load  cycles  from  FEM  simulation  of  Lane  1-1 


Lane  2-1 


Figure  6.26.  Vertical  deformation  versus  load  cycles  from  FEM  simulation  of  Lane  2-1 
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An  important  feature  for  a  model  for  granular  layers  in  pavements  is  determined  by  its 
ability  to  predict  the  shear  deformation  of  the  base  course  under  load  and  permanent  shear 
strain  under  repeated  load.  This  type  of  shearing  failure  is  representative  of  the  type  of 
behavior  seen  in  many  pavements  where  the  base  course  has  failed  (Ahlvin,  1991).  Figures 
6.27  through  6.30  show  the  evolution  of  shear  strain  through  five  load  cycles  in  Lane  1-1.  The 
shear  strain  or  principal  strain  difference  is  shown  as  a  color  fringe  plot  to  enable  the  area  of 
permanent  strain  to  be  seen.  The  shear  fringes  are  plotted  on  deformed  meshes  to  aid  in 
visualizing  the  results.  As  Load  Cycle  1  was  applied,  the  development  of  shear  strain  on  the 
order  of  0.5  %  (the  gray  shaded  region  underneath  the  loaded  area)  is  clearly  seen  in  the  base 
course  (the  2nd,  3rd,  and  4th  rows  of  elements  from  the  top)  shown  in  Figure  6.27.  After  the 
removal  of  load  cycle  1 ,  a  small  residual  shear  strain  was  developed  in  the  base  coruse  and  can 
be  seen  in  the  light  blue  shaded  region  underneath  the  loaded  area  in  Figure  6.28.  In  Figure 
6.29  the  shear  strain  under  the  5th  load  cycle  is  plotted.  The  magnitude  of  the  shear  strain 
(>0.6%)  is  greater  than  that  of  the  1st  load  cycle  due  to  the  accumulation  of  shear  strain  during 
each  of  the  five  load  applications.  The  region  experiencing  these  higher  strains  (the  darker 
shaded  region)  is  also  larger  than  that  for  one  load  cycle.  Figure  6.30  shows  the  residual  or 
permanent  shear  strain  in  the  blue,  green  and  pink  region,  after  the  removal  of  the  5th  load 
cycle.  The  permanent  shear  strain  not  only  increases  in  magnitude,  but  the  number  of 
elements  experiencing  plastic  deformation  also  increases  with  the  number  of  load  repetitions. 
The  maximum  residual  shear  strain  seen  after  the  5th  load  cycle  was  on  the  order  of  0.45%  . 
This  type  of  behavior  is  quite  representative  of  that  seen  in  pavements  subjected  to  repeated 
wheel  loads. 

Figures  6.31  through  6.34  show  the  evolution  of  shear  strain  through  five  load  cycles 
in  Lane  2-1.  The  magnitude  of  maximum  shear  strain  under  load  is  0.5%  under  one  load  cycle 
and  >1%  under  the  5th  load  cycle.  This  can  be  seen  in  the  gray  shaded  area  in  Figure  6.31  and 
the  gray/black  shaded  area  in  Figure  6.33  respectively.  The  magnitude  and  number  of 
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elements  experiencing  permanent  shear  strain  can  be  see  in  the  shaded  regions  in  Figures  6.32 
(after  1  load  cycle)  and  6.34  (after  the  5th  load  cycle).  The  maximum  residual  shear  strain 
after  5  load  cycles  is  on  the  order  of  .9%  for  Lane  2-1  which  is  almost  double  that  seen  in 
Lane  1-1.  This  would  agree  with  the  differences  in  the  two  sections  and  the  permanent 
deformation  measurements  and  predictions.  The  trend  seen  of  movement  of  material  from 
underneath  a  loaded  are  is  a  very  real  phenomenon  that  is  seen  in  pavements  under  all  types  of 
wheel  loading  conditions. 
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Figure  6.27.  Principal  strain  difference  in  Lane  1-1  (i.e.  shear  strain)  under  load  cycle  1 
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Figure  6.28.  Principal  strain  difference  in  Lane  1-1  (i.e.  shear  strain)  after  load  cycle  1 
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Figure  6.29.  Principal  strain  difference  in  Lane  1-1  (i.e.  shear  strain)  under  load  cycle  5 
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Figure  6.30.  Principal  strain  difference  in  Lane  1-1  (i.e.  shear  strain)  after  load  cycle  5 


Figure  6.31.  Principal  strain  difference  in  Lane  2-1  (i.e.  shear  strain)  under  load  cycle  1 
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Figure  6.32.  Principal  strain  difference  in  Lane  2-1  (i.e.  shear  strain)  after  load  cycle  1 
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Figure  6.33.  Principal  strain  difference  in  Lane  2-1  (i.e.  shear  strain)  under  load  cycle  5 
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Figure  6.34.  Principal  strain  difference  in  Lane  2-1  (i.e.  shear  strain)  after  load  cycle  5 
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Model  Sensitivity 


The  sensitivity  of  the  FEM  response  to  changes  in  the  Calibration  Parameters  for  the 
WES  MM  model  is  an  area  that  should  be  considered.  It  is  recognized  that  the  overall 
pavement  system  response  is  a  function  of  the  response  of  each  layer  of  the  systems  acting  as 
a  whole.  Changes  in  any  parameter  in  the  granular  material  model  will  effect  the  stress-strain 
response  and  yield  properties  of  the  material.,  although  some  parameters  are  less  important  for 
certain  types  of  materials.  This  was  demonstrated  in  Chapter  5  in  the  presentation  of  the 
MVIEWER  results.  The  real  question  is  not  how  do  changes  in  the  parameters  effect  the 
system  response,  but  how  do  changes  in  the  material  density,  compaction,  moisture  content, 
etc.,  effect  the  strength  and  stiffness  of  a  base  course  material.  A  study  to  adequately  define 
such  relationships  would  be  a  huge  effort  requiring  large  amounts  of  manpower,  laboratory 
testing,  and  time.  The  following  table  gives  some  limited  insight  into  the  question  of  the 
relationship  between  model  parameters  and  system  response.  The  global  strength  and  stiffness 
parameters,  <f>  and  G,  respectively  were  changed  form  the  original  calibration  and  the  analysis 
was  rerun.  The  deformation  in  Lane  1-1  at  the  top  of  the  base  course  under  the  5th  load,  and 
the  residual  deformation  at  the  top  of  the  base  course  after  removal  of  the  5th  load  are  shown 
for  5  cases  in  Table  6.6. 

The  overall  deformation  under  load  was  higher  with  the  lower  <j)  (Case  2)  and  lower 
with  the  higher  <j>  (  Case  1).  As  one  would  expect,  when  compared  with  the  original 
calibration,  the  overall  deformation  under  load  was  lower  in  Case  3  (higher  shear  modulus) 
and  higher  in  Case  4  (lower  shear  modulus).  The  permanent  deformations  did  not  exhibit  as 
much  sensitivity  to  changes  in  G  as  did  the  total  deformation  under  load,  which  contains  a 
large  elastic  component.  The  changes  in  response  due  to  changes  in  G  are  small  when 
compared  to  the  changes  in  <j>,  especially  in  terms  of  accumulated  deformation.  The  friction 
angle,  (|>,  was  varied  by  +/- 10%  to  arrive  at  these  results,  while  the  shear  modulus,  G,  was 
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changed  by  factors  of  2.0  and  0.5.  This  would  indicate  that  the  definition  of  the  yield  surface 
parameters  is  crucial  in  properly  modeling  the  response  of  granular  layers  in  flexible  pavement 
systems. 

The  WES  MM  model  performs  well  with  proper  calibration  and  attention  to  detail  in 
FEM  model  definition.  The  development  and  application  of  user  defined  material  models  is  a 
complex  task  that  requires  the  user  to  work  “blind”  with  an  analysis  code  that  does  not  offer 
access  to  computational  source  code.  Many  of  the  difficulties  experienced  in  applying  the 
WES  MM  model  may  be  eliminated  if  the  code  is  actually  incorporated  into  built-in  material 
libraries  instead  of  existing  as  a  user  defined  subroutine.  Such  an  endeavor  would  require  a 
collaborative  effort  with  a  commercial  analysis  code  producer  like  HKS/ABAQUS. 


Table  6.6.  Predicted  Base  Course  Deformation  in  Lane  1-1  with  Changes  in  <()  and  G. 


Friction 

Angle  <f> , 

Degrees 

Shear  Modulus. 

1  G,  ksi  (MPa) 

Predicted  Deformation  at  Top  of  Base 

Course,  mils  (mm) 

After  5m  Load 

Original  Cal 

48.0 

30  (206.8) 

165  (4.19) 

33  (0.84) 

Case  1 

52.6 

30  (206.8) 

166  (4.22) 

25  (1.00) 

Case  2 

43.2 

30  (206.8) 

188  (4.77) 

46  (1.20) 

Case  3 

48.0 

60  (413.7) 

160  (4.06) 

30  (0.76) 

Case  4 

48.0 

15  (103.4) 

193  (4.90) 

36  (0.91) 

Field  Data 

173  (4.39) 

40  (1.02) 
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CHAPTER  7:  CONCLUSIONS  AND  RECOMMENDATIONS 


CONCLUSIONS 

The  response  of  a  flexible  pavement  system  subjected  to  aircraft  loads  is  complex,  and 
accurately  predicting  the  response  of  such  a  system  requires  significant  computational 
capabilities.  For  flexible  pavements,  the  geometric  modeling  aspects  of  the  problem  are  quite 
simple,  but  the  materials  exhibit  very  complex  behavior.  The  application  of  computational 
models  in  pavements  analysis  requires  the  solution  of  many  problems  with  both  material 
constitutive  models  and  system  models.  The  materials  exhibit  viscous,  viscoelastic,  and 
plastic  response  to  loads.  Many  times  the  deformations  or  deformation  rates  are  non-linear 
functions  of  the  stress  state.  The  materials  are  often  heterogenous,  anisotropic,  and  particulate 
in  nature.  The  pavement  system  is  also  quite  complex  and  difficult  to  model.  Pavement  loads 
are  often  difficult  to  predict  over  time.  Spatial  variability  of  materials,  and  the  effects  of 
environment  and  aging  present  additional  difficulties  in  modeling  pavements.  The  purpose  of 
this  effort  was  to  develop  a  model  that  addresses  one  of  these  many  difficulties:  prediction  of 
response  of  granular  pavement  layers.  It  must  be  recognized  that  calculating  the  response  of  a 
pavement  is  of  interest  only  because  it  allows  one  to  use  it  to  predict  pavement  performance. 
The  links  between  pavement  response  and  performance  are  not  simple.  One  must  have  the 
proper  tools  to  understand  what  happens  in  pavements,  and  theoretical  models  are  needed 

The  development  of  the  WES  Multimechanical  Model  represents  a  significant 
advancement  in  the  state-of-the-art  of  flexible  pavement  material  response  modeling.  The 
essential  features  of  pavement  material  response  that  are  provided  with  the  WES 
Multimechanical  Model  include:  (1)  non-linear  elastic  response,  (2)  permanent  or  plastic 
deformations  after  yield,  (3)  cyclic  loading,  (4)  strain  softening/hardening,  and  (5)  shear 
dilatancy.  A  model  of  this  type  has  the  added  benefit  of  calibration  parameters  that  are 
physically  significant.  In  effect,  they  are  related  directly  to  the  properties  of  the  material 
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determined  from  laboratory  test  data.  The  inclusion  of  the  WES  Multimechanical  Model  for 
granular  materials  in  a  new-generation  analysis  and  design  procedure  should  provide  the 
pavements  community  with  a  tool  for  predicting  the  permanent  deformation  of  unbound  layers 
in  flexible  pavements.  The  following  conclusions  can  be  drawn  from  this  study: 

•  The  ability  to  predict  permanent  deformation  under  a  relatively  small  number 
of  load  repetitions  with  relatively  close  agreement  to  field  measurements  has  been 
demonstrated.  When  one  considers  the  long-term  effects  of  repeated  loads,  the 
analyst  must  consider  many  additional  aspects  of  the  total  pavement  system.  The 
nonlinear  plastic  response  of  the  surface  layers  and  the  subgrade  layers  are  critical  in 
understanding  the  behavior  of  a  pavement  over  time.  Variability  in  the  material 
properties  within  a  pavement  structure  is  a  systems  level  problem  that  must  be 
addressed.  The  mechanical  properties  of  granular  base  courses  can  vary  widely  from 
one  location  to  another.  The  WES  MM  model  is  quite  sensitive  to  some  of  these 
properties  and  can  result  in  differences  in  predicted  response. 

•  Historically  no  universally  accepted  rational  and  consistent  constitutive  model 
has  been  used  in  modeling  granular  pavement  materials.  A  constitutive  model  that 
can  capture  the  essential  behavior  of  pavement  materials  under  service  environments 
has  many  requirements  including  simplicity  of  calibration  and  operation,  physical 
significance  of  the  model  parameters,  and  the  ability  to  be  readily  incorporated  into 
analysis  codes. 

•  The  WES  MM  model  performs  well  in  modeling  granular  pavement  materials 
with  proper  calibration  and  attention  to  detail  in  FEM  model  definition.  It  is  essential 
that  high  quality  laboratory  test  data  be  used  to  calibrate  this  model  for  any  potential 
application.  The  selection  of  mesh  definition  for  a  pavement  section  is  more  of  a 
learned  skill  than  an  exact  science. 
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•  This  type  of  constitutive  model,  although  simple  in  formulation,  is  quite 
sophisticated  in  operation.  It  remains  a  very  high-end  analysis  tool,  which  can  be  very 
complex  in  its  application  to  pavement  analysis. 

•  The  inclusion  of  granular  base  response  models  is  critical  when  predictions  of 
pavement  behavior  under  repeated  loads  are  required.  Current  aircraft  pavement 
design  procedures  do  not  account  for  the  performance  of  the  granular  base. 

Designers  will  need  to  incorporate  criteria  for  base  course  response  in  future 
generation  analysis,  design,  and  performance  models. 

•  Older  classical  soil  models,  like  Drucker-  Prager,  appear  to  lack  the 
sophistication  required  to  properly  model  granular  pavement  materials.  The  Drucker- 
Prager  model  is  not  intended  to  handle  cyclic  loads.  The  Drucker-Prager  models  can 
not  capture  the  non-linear  pre-yield  behavior  of  granular  materials.  The  Drucker- 
Prager  model  does  not  have  the  ability  to  model  shear  dilatancy  in  materials  which 
have  been  highly  compacted. 

RECOMMENDATIONS 

The  following  recommendations  are  made  as  a  result  of  this  study: 

•  Future  generation  pavement  analysis  and  design  procedures  for  aircraft  loads 
should  include  models  capable  of  predicting  permanent  deformation  under  repeated 
loads.  The  current  pavement  design  procedures  are  capable  of  providing  reasonable 
layer  thickness  designs  for  a  wide  range  of  aircraft.  However,  when  the  task  is 
changed  to  predicting  the  performance  of  a  pavement  under  non-standard  conditions 
or  designing  a  pavement  with  non-standard  materials  then  the  older  procedures  lack 
the  sophistication  required  to  handle  that  kind  of  requirement.  This  model  would 
allow  the  designer  or  analyst  the  option  of  including  permanent  deformation  under 
repeated  load  as  criteria  in  layer  thickness  design. 


147 


•  A  database  of  test  results  necessary  for  WES  MM  model  calibration  for 
unbound  pavement  materials  should  be  developed  from  laboratory  tests  to  aid 
pavement  analysts  in  predicting  performance  of  pavements  under  repeated  loads. 
Such  a  database  would  provide  analysts  with  an  advantage  in  obtaining  values  of 
parameters  when  new  materials  are  encountered.  The  process  of  assembling  such  a 
database  would  provide  the  information  required  to  characterize  the  sensitivity  of  the 
model  and  calibration  parameters  to  changes  in  the  physical  and  mechanical  material 
properties  of  the  material. 

•  The  ABAQUS  User  Defined  Material  Model  (WES  MM)  should  be  included 
in  the  standard  material  library  for  ABAQUS  or  a  similar  code. 

•  The  WES  MM  model  should  be  the  basis  for  future  model  development  to 
include  features  such  as  partially  saturated  soils  and  time  dependent  components  for 
modeling  asphalt  cement  concrete. 

•  The  finding  of  this  study  should  be  used  in  designing  and  instrumenting  a  full 
scale  test  section  that  would  enable  the  accumulation  of  surface  and  subsurface 
permanent  deformation  of  a  pavements  under  aircraft  loads  to  be  more  accurately 
determined.  These  test  results  could  then  be  used  to  further  develop  performance 
criteria  based  on  the  WES  MM  model. 
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APPENDIX  A 

WES  MM  MODEL  UMAT  SOURCE  CODE 
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C  Last  change:  DMS  13  Jul  1999  7:57  am 

SUBROUTINE  UMAT (STRESS, STATEV, DDSDDE, SSE, SPD, SCD, 

1  RPL,  DDSDDT, DRPLDE, DRPLDT, STRAN, DSTRAN, 

2  TIME, DTIME, TEMP, DTEMP, PREDEF,  DPRED,  CMNAME, NDI , NSHR, NTENS, 

3  NSTATV, PROPS, NPROPS, COORDS, DROT, PNEWDT, CELENT, 

4  DFGRDO, DFGRD1, NOEL, NPT, LAYER,  KSPT,  KSTEP, KINC) 

INCLUDE  1 ABA_PARAM . INC  1 

CHARACTER* 8  CMNAME 

REAL* 8  STRESS (NTENS) , STATEV (NSTATV) , 

1  DDSDDE (NTENS, NTENS) , DDSDDT (NTENS ) , DRPLDE (NTENS ) , 

2  STRAN (NTENS) , DSTRAN (NTENS) , TIME (2)  ,  PREDEF (1) , DPRED (1) , 

3  PROPS (NPROPS ) , COORDS ( 3 ) , DROT (3,3),  DFGRDO (3,3), DFGRD1 (3,3) 

C!  LOGICAL  DRAINED,  Sflag(4),  Hf lag ( 4 ) , Hf lagSave ( 4 )  , 

SflagSave (4 ) 

C!  INTEGER*4  Ntimes,  iprint,  Numout,  im,  iter,  icode,  smech, 

hmech 

LOGICAL  Sf lag ( 4 ) 

LOGICAL  Hf lag ( 4 ) 

LOGICAL  Sflagd 

INTEGER*4  r 

C!  Index  for  mechanism 

INTEGER*4  counter 

REAL *8  State 
C!  Void  ratio 

REAL  *8  Qs ( 6, 4 ) ,  JQs(6,4) 

C!  Internal  shear  forces 

REAL* 8  Qh ( 4 ) ,  JQh ( 4 ) 

C!  Internal  hydrostatic  forces 

REAL*8  StateSave 
C!  Void  ratio 

REAL *8  QsSave (6,4) 

C!  Internal  shear  forces 

REAL* 8  QhSave ( 4 ) 

C!  Internal  hydrostatic  forces 

REAL *8  D (3, 3) 

Cl  Strain  Increment  tensor 

C!  REAL *8  Eps (3,3) 

Cl  Strain 

REAL* 8  Ds ( 6) 

C!  Strain  increment  vector 
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REAL* 8  Sigma (3, 3) 

C!  Stress  tensor 

REAL *8  KStress (6) 

C!  STRESS  VECTOR  USED  IN  MAIN  PROGRAM  SO  AS  NOT  TO  CONFUSE  WITH 

ABAQUS  STRESS 

REAL *8  Sigc 

C!  Confinning  stress 

REAL* 8  DeltaEps 
Strain  increment 

REAL *8  TotalEps 
Total  strain 

REAL *8  Fh,  beta,  Pe 
C!  Parameters  defining  volumetric  state 

REAL* 8  Me 

C!  Shear-volume  coupling  parameter 

REAL* 8  Cohesion 
C!  Cohesion  parameter 

REAL* 8  Gamma 

REAL* 8  C 

C!  Mohr-Coulomb  cohesion 

REAL* 8  Decay 

C!  Defines  rate  that  PhiLim  falls  with  OCR 

REAL* 8  PhiRatio 

C!  Ratio  of  maximum  and  minimum  PhiLim 

REAL* 8  PhiLim 

C!  Mohr-Coulomb  friction  angle 

REAL *8  PhiR 

C!  Friction  angle  in  radians 

REAL* 8  BulkMod 
C!  Elastic  Bulk  Modulus 

REAL* 8  ShearMod 
C!  Elastic  Shear  Modulus 

REAL *8  PhiFrac ( 4 ) 

C!  Fraction  of  PhiLim  for  each  shear  mechanism 

REAL* 8  Pfact ( 4 ) 

C!  factor  to  apportion  mean  stress  to  mechanism 

REAL* 8  Shear Ratio (4) 

C!  Shear  modulus  for  internal  mechanism 
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REAL* 8  Hlimit { 4 ) 

C!  Limit  of  internal  hydrostatic  mechanism 

REAL *8  BulkRatio ( 4 ) 

C!  Bulk  modulus  for  internal  mechanism 

REAL* 8  SPARMS ( 27 ) 


C  Counters  used  in  Do  Loops 

INTEGER*4  I,J,  IR,  IQ,  II,  IA 


C  Variables  to  calculate  Jacobian 

REAL* 8  SML_STRAIN (3,3)  ,  JAC (NTENS) 

REAL *8  JState,  BSTATE,  JSIGMA(3,3),  JSTRESS(6),  JDs(6) 
REAL *8  JACO (NTENS, NTENS) ,  ASIGMA(3,3),  blendl,  blend2 
PARAMETER (ALPHA  =  -0.00001) 

OPEN  FILE  FOR  DEBUG  DATA 


c  OPEN  (14,  FILE  =  ' . /MDUMP . OUT ' ) 

C  WRITE  (14,*)  'START  EXECUTABLE  STATEMENTS',  KSTEP,  KINC, 

NPROPS 

c  CALL  FLUSH  (14) 

c  CALL  FLUSH  (6) 

c  CALL  FLUSH  (8) 


SIGMA (1, 1) 

= 

STRESS (1) 

SIGMA(2, 2) 

STRESS (2) 

SIGMA(3, 3) 

= 

STRESS (3) 

SIGMA(1,  2) 

= 

STRESS (4) 

SIGMA(2, 1) 

= 

SIGMA(1,  2) 

SIGMA (1, 3) 

= 

0 

SIGMA (3, 1) 

= 

SIGMA (1,3) 

SIGMA(2, 3) 

= 

0 

SIGMA (3,2) 

= 

SIGMA (2,3) 

WRITE  (14, 

*) 

1  SIGMA  INITIALIZED1  ,  KSTEP,  KINC 

C  CALL  FLUSH  (14) 

C  CALL  FLUSH  (6) 

C  CALL  FLUSH  (8) 


C  PROPS (ERTIES) 


beta 

= 

PROPS (1) 

Fh 

= 

PROPS (2) 

C 

PROPS (3) 

Me 

= 

PROPS (4) 

Gamma 

PROPS (5) 

PhiLim 

= 

PROPS (6) 

Decay 

= 

PROPS (7) 

PhiRatio 

= 

PROPS (8) 

BulkMod 

= 

PROPS (9) 

ShearMod 

PROPS (10) 

PhiFrac ( 1 ) 

= 

PROPS (11) 

PhiFrac (2 ) 

= 

PROPS (12) 

PhiFrac (3) 

= 

PROPS (13) 
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C 


c 

c 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PhiFrac ( 4 )  =  PROPS (14) 
Pfact(l)  =  PROPS (15) 
Pfact (2 )  =  PROPS (16) 
Pfact (3)  =  PROPS (17) 
Pfact  (4)  =  PROPS  (18) 
ShearRatio (1)  =  PROPS (19) 
ShearRatio (2 )  -  PROPS (20) 
ShearRatio (3)  =  PROPS (21) 
ShearRatio ( 4 )  =  PROPS (22) 
Hlimit(l)  -  PROPS (23) 
Hlimit (2 )  =  PROPS (24) 
Hlimit (3)  =  PROPS (25) 
Hlimit (4)  =  PROPS (26) 
BulkRatio ( 1 )  =  PROPS (27) 
BulkRatio (2)  =  PROPS (28) 
BulkRatio (3)  =  PROPS (29) 
BulkRatio (4)  -  PROPS (30) 


WRITE  (14,*)  f  PROPS  INITIALIZED ? ,  KSTEP,  KINC 


WRITE 

(14,*) 

'  beta 

i 

f 

beta 

WRITE 

(14,*)  ' 

Fh 

i 

r 

Fh 

WRITE 

(14,*) 

'  C 

i 

r 

C 

WRITE 

(14,*) 

'  Me 

i 

t 

Me 

WRITE 

(14,*) 

'  Gamma 

i 

r 

Gamma 

WRITE 

* 

» — t 

'  PhiLim 

» 

f 

PhiLim 

WRITE 

(14,*) 

'  Decay 

i 

r 

Decay 

WRITE 

(14,*) 

'  PhiRatio 

T 

r 

PhiRatio 

WRITE 

* 

i — 1 

'  BulkMod 

? 

t 

BulkMod 

WRITE 

(14,*) 

'  ShearMod 

i 

r 

ShearMod 

WRITE 

(14,*) 

'  PhiFrac (1) 

i 

t 

PhiFrac (1) 

WRITE 

* 

t — 1 

'  PhiFrac (2) 

i 

t 

PhiFrac (2) 

WRITE 

(14,*)  ' 

PhiFrac (3) 

» 

PhiFrac (3) 

WRITE 

(14,*) 

'  PhiFrac (4) 

i 

r 

PhiFrac ( 4 ) 

WRITE 

(14,*) 

'  Pfact (1) 

i 

f 

Pfact (1) 

WRITE 

(14,*) 

'Pfact (2) 

T 

r 

Pfact (2) 

WRITE 

(14,*) 

'  Pfact (3) 

i 

r 

Pfact (3) 

WRITE 

* 

H 

'  Pfact  (4) 

i 

r 

Pfact  (4) 

WRITE 

(14,*) 

'  ShearRatio (1) 

i 

r 

ShearRatio (1) 

WRITE 

(14,*) 

'  ShearRatio (2) 

i 

r 

ShearRatio (2) 

WRITE 

(14,*) 

'  ShearRatio (3) 

i 

r 

ShearRatio (3) 

WRITE 

(14,*) 

'  ShearRatio (4 ) 

f 

f 

ShearRatio ( 4 ) 

WRITE 

(14,*) 

'  Hlimit (1) 

i 

t 

Hlimit (1) 

WRITE 

(14,*) 

'  Hlimit (2) 

f 

r 

Hlimit (2) 

WRITE 

(14,*) 

'  Hlimit (3) 

? 

f 

Hlimit (3) 

WRITE 

(14,*) 

'  Hlimit (4) 

! 

f 

Hlimit (4 ) 

WRITE 

(14,*) 

'  BulkRatio (1) 

I 

r 

BulkRatio (1) 

WRITE 

* 

^P 

i — 1 

'  BulkRatio (2) 

t 

BulkRatio (2) 

WRITE 

(14,*) 

'  BulkRatio (3) 

t 

r 

BulkRatio (3) 

WRITE 

(14,*) 

'  BulkRatio (4) 

? 

r 

BulkRatio (4 ) 

WRITE  (14,*)  1  PROPS  INITIALIZED 1 ,  KSTEP,  KINC 

WRITE  (14,*)  'COHESION  VALS ' ,  PhiLim,  Cohesion 
WRITE  (14,*)  'PhiLim  =  ',  PhiLim 
WRITE  (14,*)  'Cohesion  =  ',  Cohesion 
WRITE  (14,*)  'C  =  ',  C 
CALL  FLUSH  (14) 

CALL  FLUSH  (6) 
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CALL  FLUSH  (8) 


C  Convert  cohesion  to  a  hydrostatic  offset 

PhiR  =  PhiLim  *  3.141592/180. 


Cohesion  =  C  *  (  3.  -  SIN(PhiR)  )  *  COS (PhiR)  /  SIN(PhiR) 


C  Define  Parms 

SParms (  1)  =  beta 

SParms (  2)  =  Fh 

SParms (  3)  =  Cohesion 

SParms (  4)  =  PhiFrac(l)  *  PhiLim 

SParms (  5)  =  PhiFrac(2)  *  PhiLim 

SParms (  6)  =  PhiFrac(3)  *  PhiLim 

SParms (  7)  =  PhiFrac(4)  *  PhiLim 

Sparms (  8)  =  ShearRatio ( 1 )  *  ShearMod 
SParms (  9)  =  ShearRatio ( 2 )  *  ShearMod 
SParms (10)  -  ShearRatio (3)  *  ShearMod 
Sparms (11)  =  ShearRatio ( 4 )  *  ShearMod 
SParms (12)  =  Pfact(l) 

SParms (13)  =  Pfact(2) 

SParms (14)  =  Pfact(3) 

SParms (15)  =  Pfact(4) 

SParms (16)  =  Hlimit(l) 

SParms (17)  =  Hlimit(2) 

SParms (18)  -  Hlimit(3) 

SParms(19)  =  Hlimit(4) 

Sparms (20)  =  BulkRatio(l)  *  BulkMod 

SParms (21)  =  BulkRatio(2)  *  BulkMod 

SParms (22)  =  BulkRatio(3)  *  BulkMod 

Sparms(23)  =  BulkRatio(4)  *  BulkMod 

Sparms (24)  =  Me 

Sparms (25)  =  Decay 

Sparms (26)  =  PhiRatio 
Sparms (27)  =  Gamma 


STATE  =  STATEV (29) 

C  Load  Strain  increment  from  DSTRAN(6)  array 

0(1,1)  =  DSTRAN ( 1 ) 

D  (2 , 2  )  =  DSTRAN ( 2 ) 

D  (3, 3)  =  DSTRAN (3) 

D  ( 1 , 2 )  -  DSTRAN (4) *0.5 
D  ( 2 , 1 )  =  D(l,2) 

D  ( 1 , 3 )  =  0 
D  (3 ,  1 )  =  D  ( 1,  3) 

D (2, 3)  -  0 
D(3,2)  =  D  (2 , 3 ) 

c  Load  strain  increment  vector  Ds  from  Strain  increment  Array  D 
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Ds  ( 1 ) 

=  D ( 1, 1 ) 

Ds  (2 ) 

=  D (2,  2) 

Ds  (3) 

=  D  (3,  3) 

Ds  (4) 

=  D  ( 1 , 2 ) 

Ds  (5) 

=  D (1, 3) 

Ds  (6) 

=  D  ( 2 , 3 ) 

C  Stuff  STATE  into  BSTATE 

BSTATE  =  STATE 

Bring  Qs,  Qh,  and  Void  Ration  in  from  Statev  array. 

COUNTER=0 
DO  i-1, 6 
DO  r=l ,  4 

COUNTER-COUNTER* 1 
Qs(i,r)=  statev (COUNTER) 

END  DO 
END  Do 
DO  r*l,4 

COUNTER-COUNTER* 1 
Qh (r) =  statev (COUNTER) 

END  DO 

COUNTER-COUNTER* 1 
STATE=STATEV (COUNTER) 

c  WRITE  (14,*)  "Before”,  state,  counter  ,  ”(sv-29)", 

statev (29) 
c  DO  r=l , 4 

c  WRITE  (14,*)  ”Qh(”,R, ”) ”,Qh(r) 

c  CALL  FLUSH  (14) 

c  .  CALL  FLUSH  (6) 

c  CALL  FLUSH  (8) 

c  END  DO 

c  SUBROSand_Driver (Ds,  State,  Qs,  Qh,  Stress,  Sparms, 

Sf lag, Hflag) 

CALL  Sand_Driver (Ds,  State,  Qs,  Qh,  KStress,  Sparms, 
Sflag, Hflag) 

c  WRITE  (14,*)  "AFTER”,  state 

c  DO  r=l, 4 

c  WRITE  (14,*)  "Qh(",R, ”) ”,Qh(r) 

c  CALL  FLUSH  (14) 

c  CALL  FLUSH  (6) 

c  CALL  FLUSH  (8) 

c  END  DO 

c  CALL  FLUSH  (14) 

CALL  FLUSH  (6) 

CALL  FLUSH  (8) 


C  Fill  Qs,  Qh  and  Void  Ratio  back  into  STATE  variables 

COUNT£R=0 
DO  1-1,6 
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DO  R=l,4 

COUNTER=COUNTER+ 1 
STATEV (COUNTER)  -  Qs(I,R) 
END  DO 
END  DO 
DO  R=1 , 4 

COUNTER=COUNTER+ 1 
STATEV (COUNTER)  =  Qh(R) 
END  DO 

COUNTER=COUNTER+ 1 
STATEV ( COUNTER) =St ate 


C  Put  KStress  into  stress  variable  from  abaqus 

DO  1=1 , NTENS 

STRESS (I)  =  KSTRESS (I) 

C  WRITE (14,*)  1  STRESS  \  I,  STRESS (I) 

END  DO 


C  End  :  boss  loop 


C  CALCULATE  THE  JACOBIAN. 

C  Make  a  copy  of  Current  Qs,  Qh  and  Void  Ratio  for  a  dummy  call. 

DO  i=l, 6 
DO  r=l , 4 

QsSave(i,r)=  Qs(i,r) 

QhSave(r)  =  Qh(r) 

END  DO 
END  DO 

JState=State 

C  Zero  out  JACO (NTENS, NTENS ) 

DO  I  *  1, NTENS 
DO  J  =  1,  NTENS 
JACO (I, J)  =0.0 
END  DO 
END  DO 


C  Loop :  JACK 

DO  I  =  1, NTENS 

C  Reset  Strain  Increment  to  0.0 


DO  II  =  1,3 
DO  IQ  =  1,3 

SML_STRAIN (II, IQ)  =0.0 
END  DO 
END  DO 
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ASIGMA (1,1)  =  KSTRESS ( 1 ) 

AS I GMA (2,2)  =  KSTRESS (2) 

AS I GMA (3,3)  =  KSTRESS (3) 

ASIGMA (1,2)  =  KSTRESS (4) 

ASIGMA (1,3)  =  KSTRESS (5) 

ASIGMA (2, 3)  =  KSTRESS (6) 

ASIGMA (2,1)  =  ASIGMA (1,2) 

ASIGMA (3,1)  =  ASIGMA (1,3) 

ASIGMA (3,2)  =  ASIGMA (2,3) 

C  Set  SML_STRAIN  inc  for  partial 

IF  (I.EQ.l)  SML_STRAIN (1,1)  =  ALPHA 

IF  (I.EQ.2)  SML_STRAIN (2,2)  =  ALPHA 

IF  (I.EQ.3)  SML_STRAIN (3,3)  =  ALPHA 

IF  (I.EQ.4)  SML_STRAIN (1,2)  =  ALPHA*0.5 

IF  (I.EQ.5)  SML_STRAIN (1,3)  =  ALPHA*0.5 

IF  (I.EQ.6)  SML_STRAIN (2,3)  =  ALPHA*0.5 

SML_STRAIN (2,1)  =  SML_STRAIN (1,2) 

SML_STRAIN (3,1)  =  SML_STRAIN (1,3) 

SML_STRAIN (3,2)  =  SML_STRAIN (2,3) 

c  Load  strain  increment  vector  Ds  from  Strain  increment  Array  D 

JDs(l)  =  SML_STRAIN(1, 1) 

JDs (2)  =  SML_STRAIN (2,2) 

JDs ( 3 )  =  SML_STRAIN (3,3) 

JDs (4)  =  SML_STRAIN (1,2) 

JDs (5)  =  SML_STRAIN(1,3) 

JDs (6)  =  SML_S TRAIN (2,3) 

C  Load  the  Original  Q  Values  in  to  the  Call  Arrays  Qs,  Qh 

DO  ia=l,6 
DO  r=l, 4 

JQs(ia,r)=  QsSave(ia,r) 

JQh(r)  =  QhSave(r) 

END  DO 
END  DO 


Dummy  call  to  the  sand_driver  for  calculation  of 
Jacobian. 

CALL  Sand_Driver ( JDs,  JState,  JQs,  JQh,  JStress, 
Sparms,  Sflag,  Hflag) 


C  Map  JStress  from  Sand  Driver  into  JSigma  for  JACO 

JSIGMA (1,1)  =  JSTRESS (1) 

JSIGMA (2, 2)  =  JSTRESS (2) 

JSIGMA(3, 3)  =  JSTRESS (3) 

JSIGMA (1,2)  =  JSTRESS (4) 

JSIGMA{1, 3)  =  JSTRESS (5) 

JSIGMA (2, 3)  =  JSTRESS (6) 
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JSIGMA (2,1) 
JSIGMA  (3, 1 ) 
JSIGMA (3,2) 


JSIGMA (1,2) 
JSIGMA (1,3) 
JSIGMA (2,3) 


C  Compute  a  Jacobian  term 

IF  (I.LE.3)  THEN 
DO  J  =  1,3 

JACO ( J, I )  =  ( JSIGMA (J, I ) -ASIGMA ( J, I ) ) /ALPHA 
END  DO 

ELSE 

blendl  =  JSIGMA (1, 2) -ASIGMA (1, 2) 
blend2  =  JSIGMA(2, 1) -ASIGMA{2, 1) 

JACO (I, I)  =  (blendl+blend2 ) / (1*ALPHA) 

END  IF 
END  DO 

C  End  of  JACK  loop 

C  Zero  out  the  Jacobian  Matrix. 

DO  1=1 , NTENS 
DO  J=l, NTENS 

DDSDDE ( I , J)  =  0.0 
END  DO 
END  DO 

C  Fill  up  the  Jacobian  Matrix 

DO  I  =  1, NTENS 
DO  J  =  1,  NTENS 

DDSDDE (I, J)  =  JACO (I, J) 

END  DO 
END  DO 

END 


C  End  Program  Main 

C 

C 

******** 

c 

******** 

c  ^  ^  ^  ^  ^  ^  ^  ^  ^ . 

******** 
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SUBROUTINE 

Sand_Driver (Ds, State, Qs, Qh, Stress, Sparms, Sflag,  Hflag) 

INCLUDE  1 ABA_PARAM .  INC ' 

IMPLICIT  NONE 
C!  YIELD  SURFACE 

REAL *8  Fy 

C!  REAL *8  A (6) , B ( 6)  ,  TDOT 

REAL* 8  TDOT 

LOGICAL  Sflag (4) 

LOGICAL  Hflag (4) 

LOGICAL  Sflagd 

INTEGER* 4  r 

C!  Index  for  mechanism 

INTEGER*4  i 

C!  index  for  stress  component 

REAL* 8  Sparms (40) 

C!  Parameters 

REAL* 8  Ds ( 6) 

C!  Strain  increment 

REAL *8  DsO (6) 

C!  Null  strain  increment 

REAL *8  State 
C!  Void  ratio 

REAL* 8  Qs(6,4) 

C!  Internal  shear  forces 

REAL* 8  Qh ( 4 ) 

C!  Internal  hydrostatic  forces 

REAL* 8  Sigma 
Cl  Mean  stress 

REAL *8  S (6) 

C!  Shear  stress 

REAL* 8  SO (6) 

C!  Initial  shear  stress 

REAL* 8  Stress (6) 

C!  Stress 

REAL* 8  Sig,  SigO 
C!  Mean  stress  parameter 

REAL* 8  ShearMod 
C!  Elastic  Shear  Modulus 
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REAL *8  Fh,  beta,  Pe 

C!  Parameters  defining  volumetric  state 

REAL *8  Me 

C!  Shear-volume  coupling  parameter 

REAL *8  Cohesion 
C!  Cohesion  parameter 

REAL *8  Decay 

C!  Defines  rate  that  PhiLim  falls  with  OCR 

REAL *8  PhiRatio 

C!  Ratio  of  maximum  and  minimum  PhiLim 

REAL *8  PhiR 

C!  Friction  angle  in  radians 

REAL *8  SinPhi 

C!  Sine  of  friction  angle 

REAL* 8  Gamma 

C!  What  is  Gamma  ???????? 

REAL* 8  Phi (4) 

C!  Fraction  of  PhiLim  for  each  shear  mechanism 

REAL* 8  Ylimit ( 4 ) 

C!  Limit  of  internal  shear  mechanism 

REAL* 8  Shear (4) 

C!  Shear  modulus  for  internal  mechanism 

REAL* 8  Hlimit ( 4 ) 

Cl  Limit  of  internal  hydrostatic  mechanism 

REAL* 8  Pf act ( 4 ) 

C!  factor  to  apportion  mean  stress  to  mechanism 

REAL* 8  Bulk (4) 

C!  Bulk  modulus  for  internal  mechanism 

REAL *8  desp (6) 

C!  Plastic  shear  strain  returned  for  rth  mechanism 

REAL* 8  despt (6) 

C!  Total  plastic  shear  strain 

REAL *8  depd 

C!  Hydrostatic  strain  due  to  shear-volume  coupling 

REAL *8  dEps 

C!  Total  hydrostatic  strain  increment 

C!  Hydrostatic  Strain  increment 

dEps  =  Ds ( 1 )  +  Ds ( 2 )  +  Ds  (3) 
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n  o 


do  I  =  1,6 

desp(i)  =  0. 

END  do 

C!  Account  for  void  ratio 

State  =  (1.  +  State) *  EXP(dEps)  -  1. 
beta  =  SParms(l) 

Fh  -  SParms(2) 

Pe  -  10.** ((Fh  -  State  )*  beta) 

c  WRITE  (14,*)  "Inside",  state 

c  CALL  FLUSH  (14) 

C  CALL  FLUSH  (6) 

c  CALL  FLUSH  (8) 


C!  Fill  in  parameters 


Cohesion 

= 

SParms (  3) 

Phi (1) 

= 

SParms (  4) 

Phi  (2) 

= 

SParms (  5) 

Phi (3) 

= 

SParms (  6) 

Phi (4) 

SParms (  7) 

Shear (1) 

= 

Sparms (  8 ) 

Shear  (2 ) 

= 

SParms (  9) 

Shear (3) 

= 

SParms (10) 

Shear (4 ) 

= 

Sparms (11) 

Pfact  (1) 

= 

SParms (12) 

Pfact (2) 

* 

SParms (13) 

Pfact (3) 

- 

SParms (14) 

Pfact  ( 4 ) 

= 

SParms ( 15) 

Hlimit (1) 

= 

-  SParms (16) 

*  Pe 

Hlimit (2) 

= 

-  SParms (17) 

*  Pe 

Hlimit (3) 

= 

-  SParms (18) 

*  Pe 

Hlimit ( 4 ) 

= 

-  SParms (19) 

*  Pe 

Bulk ( 1 ) 

Sparms (20) 

Bulk (2) 

— 

SParms (21) 

Bulk (3) 

- 

SParms (22) 

Bulk (4) 

= 

Sparms (23) 

Me 

= 

Sparms (24) 

Decay 

= 

Sparms (25) 

PhiRatio 

= 

Sparms (26) 

Gamma 

= 

Sparms (27) 

ShearMod  =  Shear (1)  +  Shear (2)  +  Shear (3)  +  Shear (4) 

C!  Hydrostatic  stress  parameter 

Sig  =  Qh(l)  +  Qh (2 }  +  Qh (3)  +  Qh(4)  -  Cohesion 

Convert  friction  angle  to  yield  limit  by  building  a  principal 
stress  state  at  the  limit  and  computing  Fy  for  that  state. 

DO  r=l, 4 

PhiR= (3 . 14 15 92* Phi (r) /180. ) * 

1 (PhiRatio+ ( 1 . 0-PhiRatio) *EXP ( Decay*Sig/Pe ) ) 

SinPhi  =  SIN(PhiR) 

Stress  (-1)  =  (1.  +  SinPhi)  /  (1. -SinPhi) 

Stress(2)  =  1.0 
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Stress (3)  =  1.0 
Stress ( 4 )  =  0.0 
Stress (5)  =  0.0 
Stress (6)  =0.0 
Ylimit(r)  =  Fy (Stress) 
END  DO 

C!  Initialize  stress 

do  i=l ,  6 

S (i)  =0.0 

SO (i)  =0.0 

despt (i)  =  0.0 
desp(i)  =  0.0 
END  DO 

Sigma  =  0.0 
SigO  -  Sig 


C!  Update  each  sand  shearing  mechanism  and  shear  accumulate  stress 
DO  r  =  1,4 

C!  Save  initial  shear  stress  for  stress  dilatancy  computation 
DO  i=l, 6 

SO (i)  =  SO  (i)  +  Qs(i,r) 

END  DO 

CALL  Ammos(Ds,  Qs{l,r),  Sig*Pf act (r ) ,  desp, 

.  Ylimit(r),  Shear(r),  Sflag(r)  ) 


DO  i=l, 6 

despt (i)  =  despt (i)  +  Shear (r)  *  desp(i)  /  ShearMod 
S  (i)  =  S  (i)  +  Qs  (i,  r) 

END  DO 

END  DO 

C!  Shear  coupling  strain.  Dilation  is  positive. 

depd  =Gamma* (TDOT (SO,  despt )/( -SigO ) - 
Mc*SQRT (TDOT (despt , despt ) ) ) 


C!  Update  each  hydrostatic  mechanism  and  accumulate  hydrostatic 
stress 

DO  r  =  1,4 

CALL  Hydros (dEps-depd,  Qh(r),  Bulk(r),  Hlimit(r),  Hflag(r)) 
Sigma  =  Sigma  +  Qh(r) 

END  DO 

C!  Rescale  shear  stress  to  account  for  reduction  in  mean  stress 
Sig  =  Sigma  -  Cohesion 

DO  1=1, 6 

DsO  (I)  =0. 
s (I)  =  0. 

END  DO 


166 


DO  r  =  1,4 

DO  1=1,6 

desp(I)  =0.0 
END  DO 

CALL  Aminos  (DsO,  Qs(l,r),  Sig*Pfact  (r )  ,  desp,  Ylimit(r), 
Shear (r ) , Sflagd) 

DO  i=l , 6 

S(i)  =  S (i)  +  Qs(i,r) 

Stress (I)  =  S (i) 

END  DO 


DO  i=l, 6 


Stress ( 
END  DO 

I) 

=  S  (i) 

END  DO 

Stress (1)  = 

S 

(1) 

+  Sigma 

Stress (2 )  = 

S 

(2) 

+  Sigma 

Stress (3)  = 

s 

(3) 

+  Sigma 

c 

WRITE  (14,* 

)  "Inside  #  2",  state 

c 

CALL  FLUSH 

(14) 

c 

CALL  FLUSH 

(6) 

c 

CALL  FLUSH 

(8) 

RETURN 

END 

C  End  Subroutine  SAND_DRIVER 

C 

******** 

C 

***★***•> lr 

c 

******** 

c 


SUBROUTINE  Ammos (  Ds,  Qs,  Sig,  desp,  Ylimit,  Shear,  Sflag) 
INCLUDE  ' ABA_PARAM . INC  1 
IMPLICIT  NONE 

REAL  *8  Fy 
C!  Scalar 

REAL *8  FGrad ( 6) 

C!  Array  Giving  Gadients 
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REAL *8  TDOT 


LOGICAL  Sf lag 
REAL *8  Y,  YO 

C!  Value  of  yield  function 

REAL *8  Fc 

C!  Fraction  of  coupling  plastic  strain 

REAL *8  Sig 
C!  Mean  stress 

REAL *8  Ylimit 

C!  Limiting  value  of  yield  function 

REAL* 8  Shear 
C!  Shear  modulus 

c!  REAL *8  Gamma 

C!  Coupling  parameter 

REAL *8  Rho 

C!  Interpolation  parameter 

REAL *8  dLamda 

Cl  plastic  strain  magnitude 

REAL* 8  Ds ( 6) 

C!  Strain  magnitude 

REAL* 8  Id (6) 

C!  Identity  tensor 

REAL* 8  Qm ( 6) 

C!  Mean  stress  tensor 

REAL *8  Qs ( 6 ) ,  Qs0(6) 

C!  Shear  stress 

REAL* 8  Q (6) 

C!  Stress 

REAL* 8  des ( 6) 

C!  Shear  strain  increment  tensor 

REAL* 8  dem ( 6 ) 

C!  Volumetric  strain  increment  tensor 

REAL* 8  dQsE (6) 

C!  Elastic  strain  increment 

REAL *8  P ( 6) 

C!  Plastic  strain  direction 

REAL* 8  desp (6) 
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C!  Plastic  strain  increment  tensor 

INTEGER  I 
C!  COUNTERS 

! C  INITIALIZE  COUNTING  VARIABLES 

1=0 

C!  J=0 

C!  K=0 

C!  L=0 


C!  Identity  tensor 


Id  (1) 

=  1, 

.0 

Id  (2) 

=  1. 

.0 

Id  (3) 

=  1, 

,0 

Id  (4 ) 

=  0, 

,0 

Id(5) 

=  0. 

.0 

Id  ( 6) 

=  0. 

,0 

C!  CHECK  FOR  NON-COMPRESSION 

IF  (SIG.GE.0.0)  THEN 
SIG=-0 -  001 
END  IF 

C!  Hydrostatic  stress 


C!  BEGIN  VECTOR  COUNTER  LOOP 
DO  1=1 f  6 

Qm (I )  =  Id ( I )  *  Sig 
END  do 


C!  Save  initial  value 

Do  I  =  1,6 

QsO(I)  =  Qs ( I ) 

END  Do 

C!  Hydrostatic  increment 

DO  1=1,6 

dem(I)  =  Id ( I )  *  (Ds(l)  +  Ds(2)  +  Ds(3))/3.0 
END  DO 

C!  Shear  part 

Do  i=l, 6 

des(I)  =  Ds(I)  -  dem(I) 

END  DO 

C!  Apply  elastic  Law  with  coupling  plastic  strain 
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C!  CALL  FGradient {Q, FGrad) 

DO  1=1,  6 

dQsE(I)  =  Shear  *  des(I) 
END  DO 
DO  1=1, 6 

Qs  ( I )  =  Qs ( I )  +  dQsE(I) 
END  DO 

C!  Updated  stress 

DO  1=1,6 

Q ( I )  =  Qs (I)  +  Qm(I) 

END  DO 

C!  Trial  yield  surface 

Y  =  Fy(Q) 


c!  Adjust  stress  for  yield  condition 

IF ( Y  .GT.  Ylimit  .OR.  Y  .LE.  9.0)  THEN 
c!  Scale  back  stress 

CALL  RadialReturn (Q,  Ylimit) 

C!  Qs  =  Q  -  Qm 

DO  1=1, 6 

Qs (I)  =  Q (I)  -  Qm ( I ) 

END  DO 

c!  IF (sdump)  WRITE (13,*)  q(l),  q (2 ) ,  q(3) 

C!  Plastic  shear  strain  increment 

do  i=l, 6 

desp(i)  =(  dQsE(i)  -  (Qs (i)  -  QsO (i) )) /Shear 
end  do 

C!  Signal  that  limit  was  hit 

Sf lag  =  .True. 

ELSE 

C!  Plastic  strain  is  zero 

desp ( 1 ) =0 . 0 
desp (2) =0 . 0 
desp(3)=0.0 
desp (4 ) =0 . 0 
desp (5) =0 . 0 
desp ( 6) =0 . 0 

C!  Signal  that  limit  was  not  hit 

Sflag  =  .False. 

END  IF 

RETURN 
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o  o 


END 

C  End  of  Subroutine  AMMOS 

C! 

C!  . 

C! 

*******************************************************************'. k 

C! 

*********+*+******************************************+****+*+****+* 
Last  change:  PC  1  Apr  1999  12:49  pm 

Last  change:  PC  1  Apr  1999  12:26  pm 

C!  Subroutine  to  perform  radial  return  of  stress 

C!  point  to  yield  function  given 

C!  by  Fy(Q)  =  II  12  /I3.  A  transformation 

C!  is  first  performed  to  principal 

C!  stress  space,  then  the  return  is  performed 

C!  such  that  II  and  (Pv2-Pv3) / (Pvl-Pv3) 

C!  are  held  constant.  This  these  constraints, 

C!  Fy=Ylimit  becomes  a  cubic  equation. 

C!  The  stress  tensor  is  computed  from 

C!  the  eigen  vectors  and  adjusted  eigenvalues. 

C!  Therefore,  the  adjusted  stress  tensor  has 

C!  the  same  principal  axes,  mean  stress, 

C!  and  Lode  parameter  as  the  original  stress  tensor. 

SUBROUTINE  RadialReturn (Q,  Ylimit) 

INCLUDE  1  ABA__PARAM .  INC ' 

IMPLICIT  NONE 

LOGICAL  sdump 

LOGICAL  Reversed 

INTEGER  i,  j,  iv,  ib 
INTEGER  it 

REAL  *8  Qml,  Qm2,  Qm3 

REAL* 8  Qm 

REAL* 8  Pmag 

REAL *8  II,  12,  13,  B1 

REAL *8  A,  B,  C,  D 

REAL* 8  alpha,  beta,  gamma,  omega 

REAL *8  m (3) 

REAL* 8  f i (3) 

REAL* 8  S { 3, 3 ) 

REAL *8  Pv ( 3 ) ,  Ev (3,3) 

REAL *8  Q (6) 

REAL* 8  Ylimit,  Rmax 

C!  Initially  principal  values  not  reversed  in  order 

Reversed  =  .False. 
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C!  First  estimate  the  maximum  eigenvalue  using  Gershgorin's 

theorem 

Qml  =  Q ( 1 )  +  ABS ( Q ( 4 ) )  +  ABS(Q(5)) 

Qm2  =  ABS (Q ( 4 ) )  +Q(2)  +  ABS(Q(6)) 

Qm3  =  ABS (Q ( 5 ) )  +  ABS(Q(6))  +  Q (3 ) 

Qm  =  MAX  (Qml,  Qm2,  Qm3) 

C!  ....Compute  principal  values 

C!  Invarients  II,  12,  13 

11  =  Q  (1)  +  Q  (2)  +  Q  (3) 

12  =  Q(1)*Q(2)+Q(1)*Q(3)+Q(2)*Q(3)  -  (Q(4) **2+Q(5) **2+Q(6) **2) 

13  =  Q ( 1 ) *Q (2 ) *Q (3)  - 

Q ( 1) *Q ( 6) **2  -  Q (2 ) *Q (5 ) **2  -  Q (3) *Q (4 ) **2  + 

2 . 0*Q ( 4 ) *Q (5) *Q ( 6) 

C!  Use  Newton  iteration  to  get  largest  eigenvalue 

it  =  0 

DO  WHILE (ABS (Qm* (Qm* (Il-Qm) -12) +13) .GT.  IE-7 .AND.  it  .LE.  50) 
it  =  it+1 

Qm  =  (Qm*Qm* (2 . *Qm-Il )  +  13) / (Qm* (3. *Qm-2. *11)  +  12) 

END  DO 

C!  Compute  other  two  values  using  quadratic  obtained  from 

synthetic  division 
C!  A  =  -1.0 

B  =  II  -  Qm 
C  =  Qm  *  B  -  12 

D  =  B*B  +  4.0  *  C 

C!  D  can  be  <0  because  of  roundoff  if  there  are  repeated  roots. 

IF (D  .GT.  0.)  THEN 
D  =  SQRT(D) 

ELSE 
D=0 . 0 
END  IF 

C!  Put  in  order  of  compressive  magnitude 

Pv ( 3 )  =  Qm 

Pv (2)  =  MAX(B+D,  B-D)/2.0 
Pv(l)  =  MIN (B+D,  B-D)/2.0 


Pmag  =  MAX (  ABS(Pv(l)),  ABS(Pv(2)),  ABS(Pv(3))  ) 

C!  Check  for  null  tensor 

IF (Pmag  .LT.  l.E-12)  GOTO  777 

C!  Check  for  near-hydrostatic  conditions. 

IF (  (  ABS ( P v ( 1 ) - Pv ( 2 ) )  ) /Pmag  .LT.  l.e-3)  THEN 
IF ( (  ABS ( Pv ( 1 ) -Pv ( 3 ) )  ) /Pmag  .LT.  l.e-3)  THEN 
C!  Tensor  is  close  to  hydrostatic. 

GOTO  777 
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END  IF 
END  IF 

C!  Save  principal  values  in  normalized  form  for  use  later 

fi (1)  =  -Pv(l) /II 
f i (2 )  -  ~Pv(2)/Il 
fi (3)  =  -Pv ( 3 ) /II 

C!  Compute  principal  directions.  Note  that  by 

C!  this  point  at  least  two  eigenvalues 

C!  have  been  determined  to  be  distinct. 

C!  Order  eigenvalues  to  insure  the  first 

C!  one  is  distinct.  Note  that  they  are 

C!  now  in  order  of  magnitude.  Thus  Pv(l)  and 

C!  Pv(3)  cannot  be  equal  because  the 

C!  hydrostatic  case  has  been  ruled  out. 

IF (ABS ( Pv ( 1 ) -Pv (2 ) )  .LT.  ABS ( Pv ( 1 ) -Pv ( 3) )  .AND. 

ABS (Pv(l) -Pv(2) )  .LT.  ABS (Pv (2 ) -Pv (3) )  )  THEN 
C!  Pv(l)  and  Pv(2)  could  be  equal.  Switch  order 

Reversed  =  .true. 

A  =  Pv ( 3) 

Pv (3)  =  Pv(l) 

Pv(l)  =  A 
END  IF 

DO  i-1, 2 

IF  (i  .EQ.  1  )  THEN 

C!  First  eigenvector.  First  eigenvalue  is  distinct, 

iv  =  1 
ELSE 

C!  Pick  eigenvector  with  the  "most  distinct"  eigenvalue. 

IF (ABS ( Pv ( 1 ) -Pv ( 2 ) )  .LT.  ABS (Pv (1) -Pv (3) )  )  THEN 
iv  =  3 
ib  =  2 

ELSE 

iv  =  2 
ib=3 
END  IF 
END  IF 

C!  Set  up  the  singular  matrix 

S(l,l)  =  Q  ( 1 )  -  Pv(iv) 

S(l,2)  =  Q (4) 

S(l,3)  =  Q(5) 

S  (2 , 1)  =  S (1, 2 ) 

S  (2,2)  =  Q(2)  -  Pv (iv) 

S  (2 , 3 )  -  Q  (6) 

S  (3, 1 )  -  S(l,3) 

S(3,2)  =  S(2,3) 

S(3,3)  =  Q (3)  -  Pv (iv) 

Pmag  =  Pmag  *  Pmag 
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C!  Pick  the  appropriate  set  of  equations  for  eigenvector 

components . 

IF (ABS (S (2, 2)  *  S  (3, 3)  -  S(2,3)  *  S(3,2))/Pmag  .GT.  l.E-5) 

THEN 

D  =  S  ( 2 , 2 )  *  S  (3, 3)  -  S  (2, 3)  *  S(3,2) 

A  =  1.0 

B  =  (-S  (2,1)  *  S (3, 3)  +  S  (3, 1)  *  S (2, 3) ) /D 

C  =  (-S  (2,2)  *  S  (3, 1 )  +  S  (2 , 1 )  *  S  (2, 3) ) /D 

ELSE  IF(ABS  (S  (1,  1)  *S  (3,  3)  -S  (1, 3)  *S  (3, 1)  )  /Pmag.GT.  l.E-5)  THEN 
D  =  S(l,l)  *  S  (3, 3)  -  S (1, 3)  *  S (3, 1) 

A  =  C-S  (1,2)  *  S  ( 1 , 1 )  +  S  (3,  2)  *  S  (1,  3)  )  /D 
B  =  1.0 

C  =  (-S (1, 1)  *  S  (3,2)  +  S (3, 1)  *  S ( 1 , 2 ) ) /D 

ELSE  IF (ABS (S(1,1)*S(2,2)-S(1,2)*S(2,1)) /Pmag.GT. l.E-5)  THEN 
D  =  S  ( 1, 1)  *  S  (2 , 2 )  -  S  (1, 2)  *  S  (2, 1) 

A  =  (-S (2,2)  *  S  (1, 3)  +  S  (2, 1)  *  S (1, 3) ) /D 

B  =  (-S  (1,1)  *  S (2,3)  +  S  (2, 1)  *  S  (1, 3) ) /D 

C  =  1.0 

ELSE 

C!  Repeated  eigenvalue.  Make  a  vector  that  is  normal  to  first 

C!  and  direction  m(i)that  is  not  colinear  to  Ev(i,l) 

IF (ABS (Ev (1,1) )  .GT.  ABS(Ev(l,2))  )  THEN 
IF (ABS (Ev (1,1) )  .GT.  ABS(EV(1,3))  )  THEN 
m  (1)  *=  Ev  (3,1) 
m (2 ) =  Ev (2,1) 
m (3) =  -Ev (1,1) 

ELSE 

m(l)=  -Ev (3,1) 
m (2 ) =  Ev (2,1) 
m  (3)  =  Ev  (1, 1) 

END  IF 
ELSE 

IF (ABS (Ev ( 1 , 2 ) )  .GT.  ABS(EV(1,3))  )  THEN 
m  (1)  =  Ev  (3, 1 ) 
m(2)=  -Ev (2,1) 
m  (3)  =  Ev  (1, 1) 

ELSE 

m(l)=  -Ev (3,1) 
m(2) =  Ev ( 2 , 1 ) 
m(3) =  Ev (1,1) 

END  IF 
END  IF 

A  =  m(2)  *  Ev (3, 1)  -  m(3)  *  Ev(2,l) 

B  =  m(3)  *  Ev (1,1)  -  m ( 1 )  *  Ev(3,l) 

C  =  m(l)  *  Ev (2, 1)  -  m(2)  *  Ev(l,l) 

END  IF 

C!  Normalize  vector 

D  =  SQRT (A* A  +  B*B  +  C*C) 

Ev ( 1 , iv)  =  A/D 
Ev (2, iv)  =  B/D 
Ev(3,iv)  =  C/D 
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END  DO 


C!  Use  cross  product  to  find  third  eigenvector 

A  =  Ev(2,l)  *  Ev(3,iv)  -  Ev(2,iv)  *  Ev(3,l) 

B  =  -Ev (1,1)  *  Ev ( 3, iv)  +  Ev(l,iv)  *  Ev(3,l) 

C  =  Ev(lfl)  *  Ev(2,iv)  -  Ev(l,iv)  *  Ev(2,l) 

C!  Normalize  vector 

D  =  SQRT (A*A  +  B*B  +  C*C) 

Ev ( 1, ib)  =  A/D 
Ev (2, ib)  =  B/D 
Ev (3, ib)  =  C/D 

C!  Adjust  eigenvalues  for  yield  condition 

C!  assuming  radial  return  in  pi  plane. 

C!  The  radial  return  requires  solution 

C!  of  the  cubic  equation  that  is  obtained  by 

C!  substitution  of  fi (1) +fi (2) +fi (3) =1 

C!  and  B1  into  the  equation  for  the  yield 

C!  function.  The  root  rendering  the 

C!  largest  negative  value  (most  compressive) 

C!  is  the  correct  root.  The  cubic 

C!  is  in  the  form  of 

C!  alpha  *  Qm**3  +  beta  *  Qm**2  +  gamma  *  Qm  +  omega  =  0 

B1  =  (fi (2) -fi (3 ) ) / (fi ( 1) -fi (3) ) 

A  =  -  ( 1 . -Bl ) / (2.-B1) 

B  -  (2*Bl-l.)/(2.-Bl) 

C  =  -1 . / (2 . -Bl) 

D  =  -(l.+Bl)/(2.-Bl) 

alpha  =  B*D*Ylimit 

beta  =  B  +  B*D  +  D  + (A*D+B*C) *Ylimit 
gamma  =»  A  +  C  +  A*D  +  B*C  +  A*C*Ylimit 
omega  =  A  *  C 

C!  Use  Newton  iteration  to  get  largest 

C!  eigenvalue.  Use  approximation  from  Mohr- 

C!  Coulomb  yield  surface  as  first  guess 

Rmax  =  0.25  *  (  (Ylimit  -  5.)  +  SQRT ( ( Ylimit-9 . 0 ) * ( Ylimit-1 . 0 ) ) 

) 

Qm  =  -Rmax/ (Rmax* (Bl  +  1 . 0 ) - (Bl-2 . 0 ) ) 
it  =  0 

DO  WHILE (ABS (Qm* (Qm* (alpha*Qm+beta) -fgamma) +omega) .GT.1E-7 
.AND.  it  .LE.  50) 

it  =  it+1 

Qm  =  (Qm*Qm  *  (2.*alpha*Qm  +  beta)  -  omega)/ 

(Qm  *  (3.*alpha*Qm  +  2.*beta)  +  gamma) 

END  DO 

C!  Revised  principal  values  that  meet  yield  condition 

fi(l)  =  Qm 
fi (2)  =  A  +  B*f i { 1 ) 
fi (3)  =  C  +  D*fi(l) 

C!  Fill  back  in  to  eigenvalues 
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IF (Reversed)  THEN 

Pv(3) 

=  "fid) 

•k 

11 

Pv(2) 

=  -fi (2) 

* 

11 

Pv(l) 

ELSE 

=  -fi (3) 

* 

11 

Pv(l) 

=  “fi (1) 

* 

11 

Pv  ( 2 ) 

=  -fi (2) 

★ 

11 

Pv  ( 3 ) 
END  IF 

=  “fi (3) 

★ 

11 

C!  Rebuild  tensor  from  its  spectral  decomposition 

DO  i=l, 3 
DO  j=i, 3 

S  (i, j )  =  Pv ( 1 ) *Ev ( i ,  1 ) *Ev ( j , 1 )  + 

Pv(2) *Ev(i, 2) *Ev ( j , 2)  +  Pv (3) *Ev (i, 3) *Ev ( j , 3) 
S(j,i)  =  S  (i,  j  ) 

END  DO 
END  DO 


C!  Put  into  vector  form. 

Q  ( 1 )  =  S(l,l) 

Q  ( 4 )  =  S  ( 1 , 2 ) 

Q  (5)  =  S  ( 1 , 3 ) 

Q  (2 )  =  S  (2 , 2 ) 

Q  (6)  =  S  ( 2 , 3 ) 

Q  (3)  =  S  (3,  3) 

777  CONTINUE 


RETURN 

END 

C! 


FUNCTION  Fy (Q) 
IMPLICIT  NONE 


REAL *8  Q ( 6) 

REAL* 8  II,  12,  13 

REAL *8  Fy 

C!  Invarients  II,  12,  13 

11  =  Q ( 1 )  +  Q (2)  +  Q(3) 

12  =  Q(1)*Q(2)+Q(1)*Q(3)+Q(2)*Q(3)-(Q(4)**2+Q(5)**2+Q(6)**2) 
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13  =  Q(1)*Q(2)*Q(3)  - 

Q ( 1) *Q ( 6) **2  -  Q (2 ) *Q  ( 4 ) **2  -  G(3)*Q(5)**2  + 
2 . 0*Q ( 4 ) *Q  ( 5 ) *Q  { 6) 


Cl  Yield  Function 

Fy  =  11*12/13 

END 

C !  =============== 


C! 

C! 

C! 

C! 

C! 

C! 

C! 


*  * HYDROS*  * 


SUBROUTINE  Hydros (dEps,  Sigma,  Bulk,  Hlimit,  Hflag) 


IMPLICIT  NONE 


LOGICAL 

REAL *8 
REAL *8 
REAL *8 
REAL* 8 
REAL *8 

C!  Stress 

dSigmaE 


Hflag 

dEps 

Bulk 

Hlimit 

dSigmaE 

Sigma 

increment 
-  Bulk  *  dEps 


C!  Elastic  stress 

Sigma  =  Sigma  +  dSigmaE 


C!  Limit  condition  (note  tension — positive  convention) 

IF (Sigma  .LT.  Hlimit)  THEN 
C!  Compression  limit 

Sigma  =  Hlimit 
Hflag  =  .True. 

ELSE  IF (Sigma  . GT .  0.)  THEN 
C!  Tension  Limit 

Sigma  =  0.0 
Hflag  -  .True. 

ELSE 

Hflag  -  .False. 

END  IF 


END 

C!*  *  T  D  0  T  *  * 
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FUNCTION  TDOT (A, B) 


IMPLICIT  NONE 
C! 

C!  Function  to  compute  scalar  product  of  two  symmetric  tensors 

C!  given  in  6  vector  format 

C! 

REAL *8  A (6) ,B(6) ,  TDOT 
C! 

TDOT  =  A ( 1 )  *  B ( 1 )  +  A (2 )  *  B (2 )  +  A(3)  *  B (3)  + 

2. DO  *  (A (4 )  *  B (4 )  +  A ( 5 )  *  B (5)  +  A(6)  * 

B  ( 6)  ) 

C!  Return  TDOT  a  scalar  quantity 

RETURN 

END 


C! 

****** 

C!  This  is  a  subroutine  TO  RETURN  THE  GRADIENT  OF  A  STRESS  VECTOR 
Q(6) 

C! 

****** 


SUBROUTINE  FGradient (Q,  FGrad) 

IMPLICIT  NONE 
INTEGER  I 
REAL *8  TDOT 
REAL *8  Q(6) 

REAL *8  Iso (6) 

REAL* 8  P ( 6) 

REAL* 8  II,  12,  13 

REAL* 8  dFdll,  dFdI2,  dFdI3 

REAL *8  dlldQ ( 6) ,  dI2dQ(6),  dI3dQ(6) 

REAL* 8  Pbar, PSUM 

REAL* 8  FGrad (6) 


C!  Mean  tensor 

DO  1=1, 6 
Iso (I) =0 . 0 
END  DO 

Iso(l)  =  1. 0/3.0 
Iso (2)  =  Iso(l) 

Iso(3)  =  Iso(l) 

C!  Invarients  II,  12,  13 

11  =  Q ( 1 )  +  Q (2)  +  Q ( 3 ) 

12  =  Q(l) *Q(2)+Q(1) *Q(3)+Q(2) *  Q ( 3 )  —  (Q(4) **2+Q(5) **2+Q(6) **2) 

13  =  Q(1)*Q(2)*Q(3)  - 

Q (1) *Q (6) **2  -  Q (2 ) *Q ( 4 ) **2  -  Q(3)*Q(5)**2  + 
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2 . 0*Q ( 4 ) *Q  (5 ) *Q  ( 6) 


dFdll  =  12/13 
dFdI2  =  11/13 
dFdI3  =  -11*12/13**2 


DO  1=1,6 


dlldQ(I) 

=  Iso (I) 

END  DO 

dI2dQ(l) 

= 

Q  (2 )  +  Q<3) 

dI2dQ (2) 

= 

Q  (1)  +  Q  (3) 

dI2dQ(3) 

- 

Q  (1)  +  Q  (2 ) 

dI2dQ ( 4 ) 

= 

-2.0  *Q(4) 

dI2dQ ( 5 ) 

= 

-2.0  *Q  (5) 

dI2dQ(6) 

= 

-2.0  *Q(6) 

dI3dQ(l) 

= 

Q (2) *Q (3)  - 

Q(6) 

i  **2 

dI3dQ (2 ) 

= 

Q (1) *Q (3)  - 

Q(4] 

i  **2 

dI3dQ(3) 

= 

Q (1) *Q  (2 )  - 

Q  (53 

i  **2 

dI3dQ ( 4 ) 

= 

-2.0  * (Q (2 ) * 

Q  (4 ) 

+  Qi 

(5) *Q{6) ) 

dI3dQ ( 5 ) 

- 

-2.0  * (Q (3) * 

Q  (5) 

+  Qi 

[4) *Q(6) ) 

dI3dQ (6) 

= 

-2.0  * (Q(l) * 

Q  ( 6) 

+  Q ( 4 ) *Q(5) ) 

DO  1=1, 6 

P  (I)  - 

dFdll  *  dlldQ(I)  + 

dFdI2  *  dI2dQ (I 

+  dFdI3  *  dI3dQ(I) 


END  DO 


PSUM  =  P ( 1) +P (2 ) +P ( 3 ) 

DO  1=1, 6 

P (I)  =  P (I)  -  (PSUM) *Iso(I) 
END  DO 


PBar  =  SQRT (  TDOT(P,P)  ) 
DO  1=1, 6 

FGrad(I)  =  P(I)/Pbar 
END  DO 


RETURN 

END 


C 

********************************************************************* 
C  This  is  a  subroutine  to  intialize  the  state  variable  array  for 
ABAQUS 

SUBROUTINE  SDVINI (STATEV, COORDS , NSTATV, NCRDS , NOEL, NPT, 

1  LAYER, KSPT) 

INCLUDE  1 ABA_PARAM . INC ? 

REAL* 8  STATEV (NSTATV) , COORDS (NCRDS) 

REAL* 8  VERT, Hpart ( 4 ) 

INTEGER  counter 

C  COORDS  are  the  coordinates  of  the  point  Zero  must  be  at  the 

top  of  the  System 

NOEL  is  element  number 
NPT  is  integration  point 

LAYER  is  for  a  composite  shell  or  layered  solid 
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KSPT  is  a  section  point  with  a  curent  layer  or  section 
Hpart  is  the  fraction  of  the  bulk  modulus  in  each  mechanism 
Hydrostatic  condition  is  set  to  3*stress 

Hpart (1)  =  0.6 
Hpart (2)  =  0.38 
Hpart (3)  =0.01 
Hpart (4)  =0.01 

VERT  =  COORDS (3) *0.0868 

COUNTER=0 
DO  1=1, 6 
DO  R=1 , 4 

COUNTER=COUNTER+ 1 
STATEV (COUNTER)  =0.0 
END  DO 
END  DO 
DO  R=1 , 4 

COUNTER=COUNTER+ 1 


STATEV (COUNTER)  =  -Hpart ( R) *  1 . 0 
END  DO 


C  Set  intial  void  ratio  as  state  dependant  variable  29 

STATEV (29)  =  0.21 

RETURN 

END 
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APPENDIX  B 

SAMPLE  ABAQUS  INPUT  FILE 
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Input  File  for  Item  1-1 

Comments 

♦HEADING 

3X3  Grid  for  Lane  2 

Begin  Model  Definition 

** 

♦NODE 

Definition  of  nodes 

1, 

0.,  0. 

Node  Number,  r  coordinate,  z 

2, 

4.,  0. 

coordinate 

3, 

8.,  0. 

4, 

12.,  0. 

27, 

12.,  -260. 

1428, 

15.9827,  -260. 

• 

• 

(Lines  deleted  for  brevity) 

• 

1593, 

0.,  -260. 

1594, 

4.,  -260. 

1595, 

8.,  -260. 

** 

** 

♦ELEMENT,  TYPE=CAX4,  ELSET=AC  SURF 

1, 

1,  5,  6,  2 

Definition  of  4-node  axysynetric 

2, 

2,  6,  7,  3 

elements 

3, 

3,  7,  8,  4 

Element  number,  nodes  defining 

4, 

4,  8,  44,  10 

element 

34, 

39,  73,  74,  40 

35, 

40,  74,  75,  41 

36, 

• 

41,  75,  76,  42 

• 

• 

(Lines  deleted  for  brevity) 

♦ELEMENT,  TYPE=CAX4,  ELSET=BASE 

37, 

5,  81,  82,  6 

38, 

6,  82,  83,  7 

• 

• 

• 

(Lines  deleted  for  brevity) 

189, 

242,  276,  277,  243 

♦ELEMENT,  TYPE=CAX4,  ELSET=SUBGRADE 

217, 

100,  339,  340,  272 

218, 

272,  340,  341,  273 

• 

• 

• 

(Lines  deleted  for  brevity) 

1404, 

1591,  1595,  1427,  1393 

** 

** 

** 

End  of  Model  Definition 

** 

** 

** 

** 
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**  ACjSURF 

Material  Definition 

** 

*  SOLID  SECTION,  ELSET=AC  SURF, 

Specifies  element  properties  for  solid 

MATERIAL=AC  ELE 

1., 

** 

elements 

♦♦  BASE 

** 

*  SOLID  SECTION,  ELSET=BASE, 

MATERIAL=BASE_ELE 

1., 

** 

*♦  SUBGRADE 

** 

*  SOLID  SECTION,  ELSET=SUBGRADE, 
MATERIAL=CH6000 

1., 

** 

**  AC_ELE 

**  Date:  19-Jul-99  Time:  15:36:49 

** 

Specifies  elastic  properties  for  asphalt 

♦MATERIAL,  NAME=AC  ELE 

layer 

** 

♦DENSITY 

0.029, 

Density  in  lbs./cubic  inch 

** 

♦ELASTIC,  TYPE=ISO 

Modulus  of  elasticity  (psi),  Poisson’s 

500000.,  0.35 

ratio 

** 

**  CH6000 

**  Date:  19-M-99  Time:  15:36:49 

** 

♦MATERIAL,  NAME=CH6000 

Specifies  elastic  properties  for  subgrade 

** 

layer 

♦DENSITY 

0.023, 

Density  in  lbs./cubic  inch 

** 

♦ELASTIC,  TYPE=ISO 

12000.,  0.35 

Modulus  of  elasticity  (psi),  Poisson’s 

** 

ratio 

** 

** 

** 

♦MATERIAL,  NAME=UMAT 

Specifies  userdefined  material  for  base 

** 

course 

♦DENSITY 

0.029, 

Density  in  lbs./cubic  inch 

♦DEPVAR 

Specifies  number  for  state-dependent 

29 

variables 

** 

♦USER  MATERIAL,  TYPE=MECHANICAL, 
CONSTANTS=30,  UNSYMM 

User  defined  material  definition 

8.685,0.70,0.25,0.72, 1 .0,48.0, 1 .8,0.50 

10000., 26000.0, 0.35, 0.42,0.82,0.88, 0.9,0.77 

0.38,0.48,0.702,0.148,0.058,0.0042,0.018,0.9 

.0,1.0,0.565,0.38,0.02,0.035 

UMAT  calibration  constants 
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Comment  statements 


**  beta 
**  Fh 
**  C 
**  Me 
**  Gamma 
**  PhiLim 
**  Decay 
**  PhiRatio 
**  BulkMod 
*♦  ShearMod 
**  PhiFrac(l) 

**  PhiFrac(2) 

**  PhiFrac(3) 

**  PhiFrac(4) 

**  Pfact(l) 

♦♦  Pfact(2) 

**  Pfact(3) 

**  Pfact(4) 

**  ShearRatio(l) 
♦♦  ShearRatio(2) 
**  ShearRatio(3) 
**  ShearRatio(4) 
**  Hlimit(l) 

**  Hlimit(2) 

**  Hlimit(3) 

**  Hlimit(4) 

**  BulkRatio(l) 
♦*  BulkRatio(2) 
**  BulkRatio(3) 


End  of  UMAT  definition 


Begin  definition  of  boundary 
conditions 


(Lines  deleted  for  brevity) 


(Lines  deleted  for  brevity) 


**  BulkRatio(4) 

♦INITIAL  CONDITIONS,  TYPE=SOLUTION,  USER 
** 

** 

** 

**  AXS_Fix 
** 

♦BOUNDARY,  OP=NEW 

1, 1„  0. 

5, 1„  0. 

81,  1„  0. 


1585,  1„  0. 

1589,  1„  0. 

** 

**  BOT  FIX 
** 

♦BOUNDARY,  OP=NEW 
1427, 2„  0. 

1428, 2„  0. 


1595, 2„  0. 
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** 

**  RHS  FIX 
** 

♦BOUNDARY,  OP=NEW 

42, 1„  0. 

76, 1„  0. 

• 

• 

• 

1460, 2„  0. 

(Lines  deleted  for  brevity) 

** 

** 

** 

End  Boundary  Condition  Definition 

**  j 

**  Step  1,  Gravity 
**  LoadCase,  Geostatic 

Begin  Definition  of  Load  Steps 

** 

♦STEP,  AMPLITUDE=RAMP,  EXTRAPOLATION=NO, 
INC= 10000,  UNSYMM=YES,  NLGEOM 

Begin  Step  1:  Gravity  Load 

Application  of  Geostatic  Gravity  Load 
** 

♦STATIC 
.01,  1. 

Specify  static  analysis 

** 

Initial  time  increment,  total  time 

** 

♦CONTROLS,  PARAMETERS=FIELD, 
FIELD=DISPLACEMENT 

Specify  solution  control  parameters 

0.03, 1.0, 1.0,  , 

♦CONTROLS,  PARAMETERS=LINE  SEARCH 

6, 

♦CONTROLS,  PARAMETERS=TEME 

INCREMENTATION 

10,15,21,50,  15,,,  15,  ,6 

MM)  )1*1 
** 

** 

** 

♦DLOAD,  OP=NEW 

AC  SURF, GRAV, -2.68, 0.0, 1.0, 0.0 

Apply  distributed  gravity  load 

BASE, GRAV, -2.68,0., 1.0, 0.0 

SUBGRADE,GRAV, -2.68, 0.0, 1 .0,0.0 
** 

** 

♦FILE  FORMAT,  ASCII 
♦NODE  PRINT,  FREQ=1 

Specify  output  options 

u, 

♦NODE  FILE,  FREQ=1 

Displacements  and  rotations 

u, 

** 

♦EL  PRINT,  POS=INTEG,  FREQ=1 

s, 

E, 

Stress  and  strain 

♦EL  FILE,  POS=INTEG,  FREQ=1 

s, 

E, 
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** 

♦PRINT,  FREQ=1 

** 

♦END  STEP 

End  of  Step  1 

** 

**  Step  2,  Load 
**  LoadCase,  Default 

Begin  Step  2:  Application  of  tire  load 

** 

♦STEP,  AMPLITUDE=RAMP,  INC=10000,  NLGEOM, 
EXTRAPOLATION=NO 

Application  of  68  psi  Tire 
** 

♦STATIC 

0.001,  1,  l.E-8,  0.08 

Specify  static  analysis 

** 

Initial  time  increment,  total  time, 

** 

minimum  time  increment,  maximum 

♦CONTROLS,  PARAMETERS=FIELD, 
FIELD=DISPLACEMENT 

time  increment 

0.075,  1.0, 1.0,  , 

♦CONTROLS,  PARAMETERS=LINE  SEARCH 

Specify  solution  control  parameters 

6, 

♦CONTROLS,  PARAMETERS=TIME 
INCREMENTATION 

10,  15,21 ,50,  15,,,  15,  ,6 

Define  element  set  for  tire  load 

** 

♦ELSET,  ELSET=TIRE,  GENERATE 

1,  6,  1 

Apply  distributed  tire  load 

** 

** 

Specify  output  options 

♦♦TIRE 

** 

Displacements  and  rotations 

♦DLOAD,  OP=MOD 

TIRE,  P4,  68. 

** 

♦NODE  PRINT,  FREQ=20 

u, 

♦NODE  FILE,  FREQ=20 

u, 

|  Stress  and  strain 

i 

** 

♦EL  PRINT,  POS=CENTR,  FREQ=20 

s, 

E, 

♦EL  FILE,  POS=CENTR,  FREQ=20 

s, 

E, 

End  of  Step  1 

** 

** 

*PRINT,  FREQ=1 

** 

*END  STEP 

** 

** 

** 

** 

Begin  Step  3:  Removal  of  tire  load 

** 
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** 

** 

** 

**  Step  3,UnLoad 

**  LoadCase,  Default  Specify  static  analysis 

*STEP,  AMPLITUDE=RAMP,  INC=10000,  NLGEOM,  Initial  time  increment,  total  time, 

EXTRAPOLATION=NO  minimum  time  increment,  maximum 

Removal  of  68  psi  Tire  time  increment 

** 

♦STATIC 

0.001,  1.,  l.E-8,  0.08 

**  Apply  distributed  tire  load 

**  TIRE 
** 

♦DLOAD,  OP=MOD  End  of  Step  3 

TIRE,  P4,  0.01 

**  Load  Steps  2  and  3  are  repeated  each 

**  time  a  load  cycle  is  added 

♦END  STEP 
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APPENDIX  C 

RESULTS  OF  TRI AXIAL  COMPRESSION  TESTS 
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Axial  Strata,  pal 


Mean  Normal  Stress,  psl 


Figure  C-l.  CTC15-1 


Crushed  Limestone 
Type  610 

Traixial  Compression  at  15  psi 


Axial  Strain,  % 


Figure  C- 2.  CTC15-1 


Crushed  Limestone 
Type  610 

Triaxiaf  Corrpression  at  15  psi 


Principal  Strain  Difference,  % 


Figure  C- 3.  CTC15-1 


Crushed  Limestone 
Type  610 

TriaxiaJ  Ccrrpression  at  15  psi 
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Figure  C- 4.  CTC15-1 
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Axial  Stress,  psl  H,  Principal  Stress 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  15  psl 


kpa 

0  50  100  150  200  250  300  350  400 


;  C-  5.  CTC15-2 


Crushed  Limestone 
Type  610 

Tralxial  Compression  at  15  psl 


Axial  Strain,  % 


Figure  C-  6.  CTC15-2 
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-16  -14  -12 


-6-6-4 
Volumetric  Strain,  % 


Figure  C-  8.  CTC15-2 
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Axial  Stress,  psi 


Figure  C- 10.  CTC15-3 
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Axial  Strain,  % 


Figure  C-  11.CTC15-3 


Crushed  Limestone 
Type  610 

Triaxial  Compression  at  15  psi 


Figure  C- 12.  CTC15-3 


Volumetric  Strain,  % 


Crushed  Limestone 
Type  610 

Trlaxlal  Compression  at  30  psi 


-2  0  2  4  6  8  10  12 

Principal  Strain  Difference,  % 


Figure  C- 15.  CTC30-1 


Crushed  Limestone 
Type  610 

Trtaxial  Compression  at  30  psi 


-9-8-7-6-5-4*3-2-10  1  2 

Volumetric  Strain,  % 


Figure  C- 16.  CTC30-1 
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Mean  Normal  Stress,  psi  Jjj  Principal  Stress  Difference,  psi 


Crushed  Limestone 
Type  610 

Triaxial  Compression  at  30  psi 


C-19.  CTC30-2 


Crushed  Limestone 
Type  610 

Triaxial  Compression  at  30  psi 


-12  -10-6-6-4 

Volumetric  Strain,  % 


Figure  C- 20.  CTC30-2 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  30  psi 


Axial  Strain,  % 


Figure  C- 22.  CTC30-3 


Crushed  Limestone 
Type  610 

Triaxial  Compression  at  30  psi 


Figure  C- 23.  CTC30-3 


Crushed  Limestone 
Type  610 

Triaxial  Compression  at  30  psi 


Figure  C- 24.  CTC30-3 
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Axial  Stress,  psl 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  50  psi 


kpa 

0  100  200  300  400  500  600  700  800  900 


Figure  C- 25.  CTC50-1 


Crushed  Limestone 
Type  S10 

Traixial  Compression  at  50  psi 


Axial  Strain,  % 

Figure  C- 26.  CTC50-1 
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Figure  C- 28.  CTC50-1 
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Axial 


Figure  C- 30.  CTC50-2  Ax.aiswm.% 
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2< 


Figure  C-  33.  CTC50-3 


Crushed  Limestone 
Type  610 

TralxJal  Compression  at  50  psi 
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Axial  Strain,  % 


Figure  C-  34.  CTC50-3 


Axial  Stress,  psl  R  Principal  Stress  Difference,  psl 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  80  psi 


kpa 
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Mean  Normal  Stress,  psi 


C-  37.  CTC80-1 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  80  psl 


0  12  3 

Axial  Strain,  % 


Figure  C- 38.  CTC80-1 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  80  psi 


Axial  Strain,  % 


Figure  C-  42.  CTC80-2 
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Crushed  Limestone  Type  610 
Triaxial  Compression  at  80  psi 


Principal  Strain  Difference,  % 

Figure  C- 47.  CTC80-3 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  80  psi 


Volumetric  Strain,  % 


Figure  C- 48.  CTC80-3 


Mean  Normal  Stress,  psl 


Crushed  Limestone  Type  610 
Triaxlal  Compression  at  50  psl 


Figure  C- 49.  CTCR50-1 


Crushed  Limestone  Type  610 
Triaxlal  Compression  at  50  psl 


Figure  C- 50.  CTCR50-1 
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Figure  C- 51.  CTCR50-2 


Crushed  Limestone 
Type  610 

Triaxial  Compression  at  50  psi 


Figure  C- 52.  CTCR50-2 
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Crushed  Limestone 
Type  $10 

Hydrostatic  Compression  to  100  psi 


Axial  Strain,  % 

Figure  D- 1.  HC100-1 
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Figure  D-  2.  HC100-1 
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Figure  D-  7.  HC 100-4 


Crushed  Limestone 
Type  610 

Hydrostatic  Compression  to  100  psi 


Figure  D- 8.  HC100-4 
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Crushed  Limestone  Type  610 
Unconfined  Compression 


Crushed  Limestone  Type  610 
Unconfined  Compression 


Figure  E- 4.  UCC1 
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Crushed  Limestone  Type  $10 
Unconfined  Compression 
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Figure  E- 5.  UCC2 


Crushed  Limestone  Type  610 
Unconftned  Compression 
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Figure  E-  6.  UCC2 
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Crushed  Limestone  Type  610 
Unconfined  Compression 


Crushed  Limestone  Type  610 
Unconfined  Compression 


Volumetric  Strain,  % 

Figure  E-  8.  UCC2 
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Figure  E-  10.  UCC3 
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Figure  E- 13.  UCC4 


Figure  E-  14.  UCC4 


Crushed  Limestone  Type  610 
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E- 17.  UCC5 


Crushed  Limestone  Type  610 
Unconfined  Compression 
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Figure  E-  18.  UCC5 
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Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  F-  3.  UXE1 


Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  F-  4  .  UXE1 
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Figure  F-6.  UXE2 
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Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  F-7.  UXE2 


Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  F-8.  UXE2 
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Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 
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Figure  F-9.  UXE3 


Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  F- 10.  UXE3 


Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  F-ll.  UXE3 


Crushed  Limestone  Type  610 
Uniaxial  Strain  Test 


Figure  F-12.  UXE3 
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APPENDIX  G 

OPERATION  OF  THE  WES 
MULTIMECHANICAL  MODEL  VIEWER 


A  stand-alone  version  of  the  WES  MM  model  called  MVTEWER  was  written  to  aid  in 
determining  those  parameters  that  require  trial-and-error  methods.  MVIEWER  provides  the 
analyst  with  a  PC  compatible  platform  to  simulate  laboratory  tests  relatively  easily.  A 
discussion  of  the  MVIEWER  program  and  its  application  is  presented  in  this  appendix. 

Thirty  material  property  calibration  parameters  are  required  for  the  WES  MM  model. 
Ten  of  these  properties  are  global  (Table  G.l)  and  the  remaining  twenty  are  associated  with 
each  of  the  four  mechanisms  (Table  G.2). 


Table  G.l  Global  Properties 


Name 

Label  in  code 

Comments 

Phi 

PHILIMIT 

friction  angle 

Cohesion 

C 

cohesion 

Bulk  Modulus 

K 

Shear  Modulus 

G 

phi  ratio 

mmm 

Hydrostatic  Intercept 

Fh 

Intercept  of  Normal 

Consolidation  Line  (NCL) 

Reciprocal  of  Cc 

BETA 

Reciprocal  of  the  slope  of  NCL 

Shear-volume  factor 

Me 

shear-volume  coupling  term 

OC  factor 

Decay 

strength  reduction  term 

dilatancy  scaling  factor 

GAMMA 

Table  G.2  Mechanism  Properties 


Name 

Label  in  code 

Comments 

Strength  factor 

PHIFRAC 

scales  friction  angle 

Mean  Stress  factor 

PFACT 

scales  mean  stress 

Shear  Stiffness  factor 

SHEARRATIO 

distributes  shear  stiffness 

Compression  limit 

HLIMIT 

absolute  compression  limit 

Volumetric  Stiffness  factor 

BULKRATIO 

distributes  volumetric  stiffness 

240 


The  stand-alone  model,  MVIEWER,  was  used  to  provide  quick  feedback  during  the 
iterative  calibration  process  for  the  WES  model.  The  MVIEWER  was  compiled  using  a  LeHey 
PC  compatible  FORTRAN  77/90  Compiler.  The  MVIEWER  program  uses  either  an  ASCII 
input  file  or  an  interactive  dialogue  window  to  input  the  material  properties  and  provide  for  an 
easy  way  of  determining  the  sensitivity  of  the  WES  MM  model  to  changes  in  these  properties. 
The  main  starting  screen  for  MVIEWER  is  shown  in  Figure  G.l. 


File 


Model  Viewer 


Look  in:  I Q  M  viewer 


Ml  abaqus.dat 

M]  ctcr50.dat 

Ml  new_cal2f.dat 

Ml  calddat 

M3  ctcr50_3.DAT 

M]  new_cal2q.dat 

SB 

Ml  call.dat 

Ml  johmite.dat 

M3  new_cal2y.dat 

*|p: 

Ml  caIZdat 

Ml  new_cal.dat 

Ml  new_ca!2z.dat 

»1s 

M3  cal3.dat 

Ml  new_caI2.dat 

M3  newdatdat 

«]S 

cal4.dat 

Ml  new_cal2b.dat 

Ml  pspe_1 5 .dafc 

*|s 

Mj  cal5.dat 

Ml  new  cai2e.dat 

M3  pspe_30_1.dat 

m 

Filename:  |*.dat 

Files  of  type:  jData  Filesf'.dat) 


a 


.Open 

Cancel 


Figure  G.  1 .  Main  starting  screen  for  MVIEWER  program 


From  this  screen,  an  ASCII  data  file  containing  the  input  data  and  30  material 
properties  can  be  selected.  The  data  is  n  the  form  shown  in  Figure  G.2.  The  first  5  entries  in 
the  data  file  are  used  to  simulate  the  conditions  of  the  conventional  triaxial  test.  In  addition  to 
the  file  retrieval  method  of  inputting  data,  the  user  can  directly  type  data  into  the  appropriate 
locations  shown  in  Figure  G.3. 
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File  Edit  View  insert  Format  Help 

0.000100  0.150000  1  D 

0.010000  0.197000 

8.685000  0.700000 

0.250000  0.720000  1.000000 

48.000000  1.800000  0.500000 

10000.000000  26000.000000 
0.350000  0.420000  0.820000  0.880000 

0.900000  0.770000  0.380000  0.480000 

0.702000  0.148000  0.058000  0.004200 

0.018000  0.900000  1.000000  1.000000 

0.565000  0.380000  0.020000  0.035000  ij 

For  Help,  press  FI  j~  iNUM  .4 

Figure  G.2.  Sample  input  data  file  for  MVIEWER  program 


Figure  G.3.  Sample  input  screen  for  MVIEWER  program 
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The  MVIEWER  program  also  allows  the  analyst  the  opportunity  to  produce  plots  of 
principal  stress  difference  versus  principal  strain  difference  and  volumetric  strain  versus 
principal  strain  difference  (Figure  G.3).  Multiple  plots  from  several  runs  may  be  viewed 
together  to  aid  the  user  in  visualizing  the  effects  of  changing  the  material  properties  on  the 
stress  strain  response  of  the  model.  The  MVIWER  plot  routine  also  allows  the  user  to  plot  of 
principal  stress  difference  versus  principal  strain  difference  from  test  results  stored  in  an  ASCII 
file  (Figure  G.4)  Strains  are  given  in  %,  while  the  units  of  stress  are  determined  by  the  system 
used  in  the  calibration  (psi  or  kPa).  For  these  plots  stress  is  given  in  psi  and  strain  in  %. 


Figure  G.4.  Stress-Strain  plots  from  MVIEWER  model  results 
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Figure  G.5.  Stress-Strain  plots  from  MVIEWER  model  results  (longer  upper  line)  and  test  data 
(lower  shorter  line) 
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APPENDIX  H 

DETERMINATION  OF  STRENGTH  PARAMETERS 
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The  Mohr  circle  of  stress  provides  a  convenient  method  of  analyzing  two-dimensional 
stress  states.  In  order  to  apply  the  method,  the  values  and  directions  of  the  principal  stress 
must  be  known.  In  the  case  of  conventional  triaxial  tests  of  soils  the  applied  stresses  are  the 
principal  stresses.  The  axial  stress  is  the  maximum  principal  stress  (cti)  and  the  confining 
stress  is  the  minimum  principal  stress  (a3).  The  maximum  shear  stress  has  the  coordinates  of 
(s,  t)  as  shown  in  Figure  H.  1 . 


Figure  H.  1 .  Mohr  circle  of  stress  for  a  conventional  triaxial  compression  test 

In  the  case  of  plastic  analysis  of  soils  behavior,  the  Mohr  circle  containing  the  normal 
and  shear  stresses  at  failure  is  a  limiting  circle.  Limiting  circles  at  different  values  of  normal 
stress  will  all  touch  at  a  common  tangent,  which  is  called  a  failure  envelope  (Figure  H.2).  The 
equation  of  this  failure  envelope  is  referred  to  as  Coulomb’s  equation: 

t  =  c  +  <j  tan  (f)  (H.l) 
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Where: 

x  =  shear  stress 

c  =  cohesion 

a  =  normal  stress 

<(>  =  angle  of  internal  friction 

A  line  drawn  through  the  point  of  maximum  shear  stress  (s,t)  for  a  series  of 
conventional  triaxial  compression  tests  will  produce  a  maximum  stress  point  failure  envelope. 
The  equation  of  this  line  is  given  as: 

t  =  a  +  s  tan  a  (H.2) 

Where: 

t  =  1/2  (ci-g3) 
s  =  1/2  (ai+a3) 
a  =  intercept  (c  cos  a=  a) 
a  =  friction  angle  (sin  (j)  =  tan  a) 
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Figure  H.2.  Failure  envelopes  from  Mohr’s  circle  of  stress  for  two  conventional  triaxial 
compression  tests 
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The  version  of  the  WES  Multimechanical  Constitutive  Model  (WES  MM)  used  in  the 
research  reported  in  the  main  body  of  this  dissertation  was  originally  formulated  for  full  three 
dimensional  (3D)  analyses.  The  model  was  simplified  to  operate  in  a  two  dimensional  axis- 
symmetric  case.  The  laboratory  and  field  tests  analyzed  were  well-suited  to  an  axisymmetric 
analysis.  In  future  analyses  the  investigation  of  multiple  wheel  response  and  moving  loads 
will  require  that  the  pavement  system  to  be  modeled  in  a  full  three  dimensional  setting. 

Since  the  original  formulation  of  the  WES  MM  was  3D,  it  was  relatively  simple  to  set 
the  model  back  to  operate  with  a  3-D  8-node  isoparametric  brick  element.  In  order  to 
demonstrate  the  effectiveness  of  the  model  in  three-dimensional  analysis,  a  single  1-in.  cubical 
element  was  subjected  to  the  same  stress  path  as  the  50-psi  conventional  triaxial  compression 
test.  The  element  was  subjected  to  the  3-D  equivalent  of  the  load  and  boundary  conditions  in 

the  axisymmetric  analysis  presented  in  Chapter  6.  The  horizontal  stresses  (C72  and  CT3)  are 

held  at  50  psi,  while  the  vertical  stress  (CT 1)  was  increased  until  a  vertical  strain  (81)  of 

approximately  5%  was  achieved.  The  element  and  the  boundary  conditions  are  shown  in 
Figure  1.1. 

The  laboratory  test  results  and  FEM  predictions  are  shown  in  Figure  1.2  and  1.3.  The 
3D  analysis  is  slightly  stiffer  at  high  strain  levels  that  the  2D  analysis.  This  can  be  attributed 
to  the  differences  in  element  formulation  in  ABAQUS  and  small  differences  in  convergence 
criteria.  The  maximum  difference  between  the  2D  and  3D  predicted  stress  is  only  3%.  The 
application  of  the  WES  MM  to  3D  problems  is  an  area  for  future  exploration  with  many 
applications  in  pavement  analysis. 
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Figure  1.1.  3D  element  under  triaxial  compression  loading 


Crushed  Limestone  Type  610 
Triaxial  Compression  at  50  psi  (344.7  kPa) 


Principal  Strain  Difference,  % 


Figure  1.2.  Laboratory  test  results  and  FEM  predictions  for  a  50  psi  triaxial  compression  test 
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Crushed  Limestone  Type  610 
Triaxial  Compression  at  30  psi  (206.8  kPa) 
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Figure  1.3.  Laboratory  test  results  and  FEM  predictions  for  a  30  psi  triaxial  compression  test 
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